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Preface 


This  thesis  examines  the  problem  of  wave  propagation  in  a  randomly  inhomogeneous 
medium  from  the  standpoint  of  Feynman  integration.  Also  known  as  path  integration  or 
functional  integration,  Feynman  integration  oftentimes  presents  the  uninitiated  reader  with 
conceptual  difficulties;  therefore,  this  thesis  was  written  with  a  tutorial  style  in  mind.  Moreover, 
there  are  numerous  appendices  and  far  greater  detail  and  verboseness  in  the  main  chapters  than 
would  normally  be  considered  necessary  in  a  thesis;  however,  it  is  hoped  that  excessive  verbosity 
will  be  an  asset  in  this  case. 

Kyle  Hunter 
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Abstract 

The  problem  of  calculating  the  scintillation  index  of  an  atmospherically  propagating 
spherical  wave  is  examined.  The  fourth  statistical  moment  of  the  complex  field  is  obtained  by 
using  Feynman  (path)  integral  techniques  applied  to  the  stochastic  parabolic  equation.  The 
general  trajectory  of  each  Feymnan  integral  is  approximated  by  a  truncated  Fourier-sine  series 
and  the  infinite-fold  integration  of  the  Feynman  integral  is  reduced  to  a  three-fold  Riemann 
integral  which  is  shown  to  match  results  derived  under  different  assumptions. 


WAVE  PROPAGATION  IN  A  RANDOMLY  INHOMOGENEOUS  MEDIUM  -  A 

STUDY  OF  THE  PROBLEM 


/.  Introduction  to  the  Problem 


The  Overview 

W  ve  propagation  a  randomly  inhomogeneous  medium,  for  example  monochromatic 
light  propagating  in  the  atmosphere  results  in  changes  in  both  the  phase  and  the  amplitude  of  the 
propagating  wave.  If  the  fluctuations  of  the  amplitude  and  phase  are  weak,  then  it  is  possible  to 
describe  the  characteristics  of  the  amplitude  and  phase  using  perturbation  techniques.  In  this 
case,  the  observed  changes  of  the  phase  of  the  wave  in  an  observation  plane,  some  distance  from 
the  wave’s  origin,  will  be  small.  The  region  over  which  an  observation  plane  may  be  placed,  and 
still  facilitate  the  measurement  of  small  phase  fluctuations,  is  known  variously  as  the  weak 
fluctuation  regime,  the  perturbational  regime,  or  the  asymptotic  regime.  Moreover,  as  the 
propagation  distance  increases  or  the  strength  of  the  random  inhomogeneities  increases,  the 
intensity  changes  in  the  observation  plane  of  the  propagating  wave  will  increase.  For  sufficiently 
strong  random  inhomogeneities  in  the  atmosphere,  the  intensity  variations  will  attain  an 
extremum  and  the  scintillation  index,  the  normalized  intensity  variance,  may  exceed  unity.  If  the 
scintillation  index  exceeds  unity,  then  the  propagation  distance  or  region  over  which  scintillation 
index  does  so  is  known  as  the  strong  fluctuation  regime.  Finally,  if  the  propagation  distance  is 
increased  further  still  or  the  random  inhomogeneities  become  stronger  yet,  then  the  intensity 
fluctuations  will  begin  to  decrease  and  asymptotically  saturate  at  a  fixed  level.  This  region  of 
asymptotic  saturation  is  known  as  the  saturation  regime.  Wave  propagation  in  random  media  are 
characterized  by  the  medium’s  statistical  characteristics.  This  implies  that  various  statistics,  and 
particular  moments  must  be  obtained.  These  can  be  found  by  averaging  over  stochastic  complex 
wave  functions  obtained  by  Feynman  (path)  integration. 
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The  Problem 


Functional  integration,  by  its  vary  nature,  requires  multiple  integrations  over  a  continuum 
of  elements  of  a  function  space.  There  is  presently  a  method  to  approximate  a  functional  integral, 
and  therefore  the  integration  continuum,  by  using  a  finite  subset  of  an  infinite  orthogonal  function 
expansion  to  approximate  the  functional  integrand.  The  resulting  approximation  can  then  be 
numerically  integrated  to  obtain  the  necessary  complex  amplitude.  The  accuracy  of  some  of  the 
initial  subsets,  a  finite  number  of  terms  from  a  Fourier  series,  looks  ver^  promising.  The 
application  of  these  techniqes  to  the  proplem  of  wave  propagation  in  a  randomly  inhomogeneous 
medium  is  a  fairly  recent  development  and  not  all  of  the  ramifications  are  at  present  fully  known 
or  understood. 

The  Research  Questions 

During  the  course  of  research  for  this  thesis,  it  desired  to  obtain  the  answers  to  some 
relevant  questions. 

1.  Using  only  one  sinusoid,  the  first  term  in  a  Fourier-sine  series,  what  accuracy  does  it 
provide  when  used  to  approximate  functional  integration? 

2.  What  is  the  computational  complexity  of  this  approximation? 

3.  Does  the  characteristics  of  this  approximation  dictate  the  use  of  a  particular  numerical 
integration  scheme? 

These  questions  are  obviously  not  all  inclusive.  In  fact,  it  will  be  necessary  to  drastically 
reduce  the  potential  research  extent  in  order  to  obtain  any  meaningful  results.  This  will  be  done 
by  making  certain  assumptions. 

The  Research  Assumptions 

As  it  will  be  explained  in  greater  depth  in  Chapter  2,  all  propagating  electromagnetic 
waves  must  satisfy  the  wave  equation  which  is  derived  from  Maxwell’s  equations.  However,  this 
derived  equation  is  a  second  order  in  space,  second  order  in  time,  partial  differential  equation.  If 
the  situation  is  further  compounded  by  requiring  propagation  in  a  randomly  inhomogeneous 
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medium,  then  the  wave  equation  becomes  a  second  order  stochastic  partial  differential  equation. 
By  making  several  prudent,  realistic  assumptions,  it  is  possible  to  simplify  the  solution  to  the 
wave  equation,  a  Feynman  integral,  into  a  solution  that  is  less  complicated.  In  order  to  make  this 
simplification,  it  will  be  necessary  to  make  the  following  assumptions: 

1.  The  parabolic  approximation  to  the  stochastic  Helmholtz  equation  holds. 

2.  The  propagating  wave  is  monochromatic. 

3.  The  wavelength  of  the  wave’s  carrier  is  in  the  optical  regime. 

4.  The  characteristic  size  of  the  random  inhomogeneities  of  the  propagation  medium  greatly 
exceeds  the  wavelength  of  the  carrier. 

5.  The  velocity  of  light  will  greatly  exceed  the  maximum  velocity  of  the  fluctuations  which 
implies  nonrelativistic  motion. 

6.  Near  field  effects  are  not  relevant  as  they  are  small  and  we  will  consider  propagation 
distances  which  exceed  the  characteristic  length  of  the  near  field. 

7.  Although  the  Feynman  integral  methods  was  borrowed  from  quantum  mechanics,  quantum 
mechanical  effects  are  assumed  to  be  inconsequential  in  all  limiting  cases  of  the  functional 
integrals. 

Assumption  1,  the  parabolic  approximation,  will  require  that  back  scattering  effects  as 
well  as  large  scattering  angles  from  the  mean  direction  of  propagation,  be  ignored.  As  it  is  shown 
in  Appendix  A,  this  necessitates  that  the  complex  amplitude  of  a  propagating  wave  vary  slowly  as 
compared  to  the  mean  wavelength  of  the  carrier. 

Assumption  2  does  not  excessively  limit  the  applicability  of  the  results  of  this  thesis;  it  is 
possible  to  decompose  a  general  propagating  wave  into  its  spatial  Fourier  components  and  then 
consider  each  spatial  frequency  separately;  however,  this  thesis  will  not  cover  such  an  analysis. 

Assumption  3  will  allow  effects  of  wave  depolarization  due  to  the  random 
inhomogeneities  to  be  ignored  because  the  depolarization  will  be  small. 

Assumption  4  is  realistic  as  the  size  of  the  inhomogeneities  in  Earth’s  atmosphere  ranges 
from  several  millimeters,  the  inner  scale,  to  about  a  hundred  meters,  the  outer  scale. 
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Assumption  5  requires  a  physically  realizable  propagation  medium  such  as  the  earth’s 
atmosphere.  Notice  this  assumption  might  not  be  valid  for  optical  wave  propagation  in  a 
turbulent  plasma. 

Assumption  6  states  that  the  near  field  effects  may  be  ignored.  The  Huygens-Fresnel 
principle  provides  a  rigorous  solution  to  the  parabolic  wave  equation  for  positive  propagation 
distances  which  exceed  the  near  field  regime.  This  will  be  covered  in  more  detail  in  Chapter  3. 

Assumption  7  requires  the  consideration  of  only  classical,  that  is,  non-quantum 
mechanical  phenomena. 

77ie  Scope 

This  thesis  will  investigate  the  problem  of  wave  propagation  in  random  media  and  some 
of  the  equations  that  will  be  generated  by  this  investigation.  Additionally,  an  accleration 
technique  which  can  be  used  to  solve  the  final  resulting  equation  will  be  provided. 

The  Research  Methodology 

The  numerical  integration  of  continuous  functions  is  a  field  of  study  in  itself.  The 
integrand  functions,  which  will  be  considered  in  this  thesis,  are  highly  oscillatory;  therefore, 
many  of  the  more  favored  numerical  integration  schemes,  such  as  Gaussian  quadrature,  will  not 
achieve  acceptable  accuracy.  Thus,  in  general,  one  finds  that  for  many  currently  available 
numerical  integration  routines  that: 

1.  The  proposed  numerical  integration  scheme  may  not  converge  within  the  prescribed 
number  of  iterations. 

2.  The  proposed  numerical  integration  scheme  may  not  converge  at  all.  Some  integration 
schemes  recursively  subdivide  every  interval  and  compare  the  results  of  each  subdivision. 
If  the  solution  series,  which  is  generated  by  comparing  the  integral  over  an  interval  with  all 
the  recursive  integrals  of  that  interval,  is  bounded,  but  fails  to  converge,  or  worse  yet, 
diverges  to  infinity,  then  the  solution  and  the  associated  method  of  solution  are  useless. 

3.  The  proposed  numerical  integration  scheme  rapidly  converges  to  an  incorrect  solution.  This 
is  possible  if  each  subdivision  matches  some  pathological  characteristic  of  the  integrand. 
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The  Research  Benefits 

As  this  thesis  is  directed  primarily  towards  a  study  of  the  problem  of  analyzing  wave 
propagation  in  a  randomly  inhomogeneous  medium,  the  benefits  of  this  research  will  primarily  be 
to  enable  the  reader  of  this  thesis  to  better  understand  the  problem  and  its  characteristics.  That  is, 
by  casting  a  physical  problem  into  an  equivalent  mathematical  form,  it  is  possible  to  obtain  a 
deeper  understanding  of  said  problem.  The  Feynman  integral  method  will  provide  an  intuitively 
pleasing  way  to  analyze  the  problem;  therefore,  it  will  allow  the  reader  to  obtain  greater  insight 
into  the  problem  of  wave  propagation  in  a  randomly  inhomogeneous  medium.  This  basic 
understanding  of  a  difficult  problem  is,  and  should  always  be,  a  goal  of  basic  theoretical  research 
into  said  problem. 

The  next  chapter  contains  a  literature  review  of  wave  propagation  analysis  using 
Feynman  integration,  as  well  as  other  currently  favored  analysis  methods. 
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//.  Literature  Review 


Justification 

Imagine  a  shinny  penny  on  the  bottom  of  a  clear,  dark-bottomed  swimming  pool  on  a 
calm  day.  Standing  on  the  edge  of  the  pool,  one  can  easily  see  that  penny.  If  you  toss  in  a  rock, 
the  penny’s  image  will  dance  about,  and  it  may  even  disappear  entirely.  The  atmosphere  is 
obviously  not  that  swimming  pool,  however  both  the  hypothetical  pool  and  the  atmosphere 
randomly  distort  optical  images  by  distorting  propagating  optical  waves.  That  is,  the  propagating 
optical  wave  is  distorted  randomly  in  both  phase  and  in  amplitude  by  the  medium  in  which  it  is 
propagating.  It  is  natural  to  inquire  as  to  characteristics  of  the  propagating  optical  wave.  It  is  the 
characteristics  of  an  atmospherically  contaminated  propagating  optical  wave  with  which  this 
thesis  will  be  concerned. 

The  atmosphere  is,  as  a  first  approximation,  a  chunk  of  air  with  varying  densities 
generated  by  thermal  heating  from  the  sun.  Large  sections  of  air  warm,  rise,  and  eventually  break 
into  smaller  sections.  The  smaller  sections  disperee  into  even  smaller  air  pockets.  These  air 
pockets  can  range  in  size  from  a  few  millimeters,  the  inner  scale,  to  several  hundred  meters,  the 
outer  scale  (17:388,390).  It  is  the  motion  of  the  air  pockets,  along  with  their  varying  relative 
indices  of  refraction,  that  cause  the  randomization  of  the  amplitudes  and  phases  of 
atmospherically  propagating  optical  waves.  Because  the  atmosphere  introduces  these 
randomizations  into  propagating  optical  waves,  the  characteristics  of  the  optical  waves,  their 
amplitude  and  phase,  are  unpredictable.  Therefore,  it  is  the  statistical  characteristics  of  the 
propagating  optical  waves  which  must  be  considered.  Moreover,  the  method  by  which  the  optical 
wave’s  statistics  are  obtained  will  be  dictated  by  the  fluctuating  statistical  characteristics  of  the 
medium,  and  the  distance  through  which  the  optical  wave  has  propagated.  There  are  several 
mathematical  methods  which  can  be  used  to  either  calculate  or  approximate  the  statistical 
characteristics  of  an  optical  wave  propagating  through  a  randomly  inhomogeneous  medium,  such 
as  the  atmosphere.  These  mathematical  methods  will,  in  turn,  have  conditions  under  which  they 
are  applicable. 

The  choice  of  the  particular  mathematical  analysis  method  is  often  complicated  by  the 
fact  that  the  characteristics  of  the  optical  wave  can  influence  the  choice  of  the  mathematical 
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method  used  to  analyze  them.  In  order  to  gain  a  more  thorough  understanding  of  the  especially 
powerful  approach  of  functional  integration,  it  is  necessary  to  be  familiar  with  some  of  the  other 
methods  which  are  currently  mentioned  in  the  literature  along  with  their  regions  of  applicability 
and  any  constraints  to  their  use.  To  do  this,  it  will  be  necessary  to  start  at  the  beginning  with  the 
wave  equation. 

In  free  space,  a  propagating  monochromatic  optical  wave  must  obey  the  wave  equation: 

(2.1) 

where!/ (r,  t)  is  the  optical  field,  a  function  of  space  and  time,  ^  ^  ^  is  the 

position  vector  (i,  y,z),c  is  the  speed  of  light,  t  is  time.  A  closed  solution  of  Eq.  2.1  is,  except 
for  trivial  or  highly  symmetric  cases,  extremely  difficult,  if  not  impossible,  to  obtain.  Therefore, 
numerical  solutions  to  Eq.  2. 1  must  be  obtained.  However,  the  solution  of  Eq.  2.1  requires  a 
global  approach.  That  is,  the  solution  must  be  obtained  in  all  volumes  of  interest,  subject  to  all 
the  boundary  conditions  simultaneously  (23:1).  It  is  possible  to  reduce  the  computational 
complexity  of  Eq.  2.1  by  considering  a  reduced  class  of  problems.  By  considering  cases  with 
negligible  backscatter  and  small  propagation  angles  about  the  mean  direction  of  propagation,  Eq 
2.1  will  be  simplified.  This  simplification,  which  results  in  a  new  equation,  is  often  called  the 
parabolic  approximation  to  the  wave  equation.  These  characteristics  of  this  simplified  equation 
will  be  explored  in  the  next  section. 

Parabolic  Equation 

The  derivation  of  the  parabolic  approximation  to  the  wave  equation  for  an  optical  wave 
propagating  in  the  z  direction  can  be  found  in  Appendix  A.  It  is  shown  in  Appendix  A  that 
Eq.  2.1  may  be  reduced  to 

(2;*^  +  vi  +  f  ip,  z)Mp,  z)  =  0  (2.2) 

where  u{p,  z)  is  the  complex  amplitude  of  the  optical  wave,  ^  i  is  \/M,  fc  is  ^ 
and  X  is  the  carrier  wavelength,  e  (p,  z)  is  the  fluctuating  component  of  the  medium’s  dielectric 
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permeability,  p  =  (i,  y),  namely  p  is  in  the  plane  transverse  to  the  mean  direction  of  propagation, 
z. 

The  simplification  of  Eq.  2.1  to  Eq.  2.2  is  a  direct  result  of  considering  the  propagation  of 
the  optical  wave  in  small  angles  about  the  mean  direction  of  propagation.  The  small  angle 
simplification  which  results  in  Eq.  2.2  has  been  known  under  various  names.  Some  of  these 
names  include:  "small  angle  approximation",  "paraxial  approximation",  "the  standard  parabolic 
approximation"  (14:838),  and  the  randomly  occurring  implied  "usual  approximation"  ,  where  the 
article’s  author  assumes  a  certain  level  of  reader  sophistication  (2:516).  The  reason  for  this 
interest  in  small  angles  about  the  mean  propagation  direction  is  that  some  random 
inhomogeneous  media,  under  certain  conditions,  exhibit  only  small  wave  field  scattering  angles. 
Therefore,  the  wave  equation  may  be  reduced  to  the  parabolic  wave  equation.  However,  this 
small  angle  restriction  is  more  fundamentally  a  necessary  condition  to  maintain  the  validity  of  the 
approximation  of  Eq.  2.1  by  Eq.  2.2.  Specifically,  the  backscatter  of  the  optical  wave,  or 
acoustical  wave  for  that  matter,  is  considered  to  be  negligible  (10:895).  The  interested  reader  will 
find  Appendix  A  has  further  details. 

The  parabolic  approximation  to  the  wave  equation  is  an  active  area  of  research  for 
numerous  authors  including  (29:212),  (10:894),  (8:171),  (7:1224),  (6:38-3),  (12:297),  and 
(13:1355).  Furthermore,  authors  in  other  fields  including  acoustics  and  quantum  mechanics  are 
also  extremely  interested  in  this  equation  because  of  the  existence  of  equations  which  are 
functionally  similar  to  Eq  2.2.  In  acoustics,  the  propagation  of  sound  in  the  ocean  is  modeled  as 
an  inhomogeneous  second  order  partial  differential  equation  (23:6).  Therefore,  the  general 
solution  to  Eq  2.2  and  the  acoustics  problem  will  be  similar.  In  quantum  mechanics,  the  motion 
of  elementary  particles  in  a  conservative  field  is  modeled  by  an  inhomogeneous  partial 
differential  equation  known  as  the  Schrodinger  wave  equation  (29:225),  (16:48).  This  equation  is 
equivalent  to  the  parabolic  approximation  to  the  wave  equation  and  therefore  general  solutions 
will  be  similar  to  both  equations.  It  is  indeed  fortunate  that  these  diverse  fields  possesses  similar 
equations,  because  the  solution  obtained  in  one  field  of  study  may  be  used  in  the  other  research 
fields,  subject  to  the  appropriate  boundary  conditions  (29:225,226),  (23:7,8,18,19).  There  is  a 
plethora  of  different  names  for  Eq.  2.2  in  the  current  literature.  Eq.  2.2  is  also  known  as  the 
parabolic  approximation  to  the  stochastic  Helmholtz  equation,  the  parabolic  approximation  to  the 
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wave  equation  (14:838),  reduced  scalar  Helmholtz  equation  (12:297),  complex  heat  equation 
(7:1224),  diffusion  equation  (16:49). 

As  it  was  stated  earlier,  Eq.  2.1  usually  doesn’t  posses  a  closed-form  solution,  and  the 
random  effects  of  the  atmosphere  necessitate  the  study  of  statistical  characteristics,  rather  than, 
the  deterministic  characteristics  of  the  optical  wave.  However,  by  considering  Eq.  2.2  and  the 
statistics  of  the  atmosphere  (more  generally  the  medium  in  which  propagation  occurs),  it  is 
possible  to  apply  different  mathematical  techniques  to  the  solution  of  Eq.  2.2.  The  next  sections 
will  consider  the  various  mathematical  methods. 

Method  of  Geometrical  Optics 

The  method  of  geometrical  optics  is  the  lowest  approximation  to  the  wave  field 
amplitude  in  a  random  inhomogeneous  medium.  In  this  multiple  scattering  method,  the  accuracy 
of  the  solution  increases  with  decreasing  wavelength  (10:894).  Furthermore,  the  range  of  validity 
for  geometric  optics  can  be  large  (25:2).  Additionally,  the  method  facilitates  the  analysis  of 
multiple  effects,  including  refraction,  and  is  the  basis  for  the  method  of  smooth  perturbations  and 
the  Markov  approximation  method  (25:2). 

Geometrical  optics,  the  method  with  which  most  individuals  are  familiar,  is  the  method 
whereby  optical  fields  are  modeled  by  straight  line  or  ray  constructions.  This  method  is  the 
simplest  of  the  analytical  methods  used  to  obtain  a  solution  to  the  problem  of  wave  propagation 
in  a  randomly  inhomogeneous  medium.  The  simplicity  is  achieved  by  ignoring  diffractional 
effects.  The  simplicity  is  attained  when  the  umbra  and  penumbra  of  shadows  are  determined  by 
the  geometry  of  transparent,  translucent,  and  opaque  objects  in  the  propagation  medium.  Because 
of  this  simplicity,  the  geometric  optics  method  can’t  predict  the  effects  of  diffraction  nor  can  it 
take  them  into  account.  Thus,  simplicity  is  traded  for  accuracy;  "The  method  is  simple  in  that, 
unlike  the  method  of  smooth  perturbations  and  Markov  approximation,  it  does  not  take  into 
consideration  diffraction  effects,  and  therefore  is  not  as  versatile  as  the  other  two  methods" 

(25:1).  However,  the  method  has  certain  advantages  over  other  methods.  Specifically, 
geometrical  optics  permits  multiple  effect  modeling  and  the  modeling  of  effects  that  are  more 
difficult  to  describe  using  other  methods.  Additionally,  certain  propagation  effects  will  continue 
to  maintain  validity  beyond  the  geometrical  optics  (validity)  regime  (25:2).  The  elimination  of 
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diffractional  effects  may  be  justified  under  certain  geometries;  however,  when  diffraction  must  be 
considered,  the  method  of  smooth  perturbations  provides  a  possibly  better  solution  framework. 


Method  of  Smooth  Perturbations 

This  single  scattering  method  is  one  of  the  older  methods  still  being  used  to  solve  Eq.  2.2. 
This  method,  also  known  as  either  Rytov’s  method,  or  the  asymptotic  method  in  the  current 
literature,  uses  first  order  perturbation  methods  on  the  logarithm  of  the  optical  wave’s  complex 
amplitude  to  compute  an  approximate  solution  to  Eq.  2.2.  This  method  may  be  used  whenever 
the  scintillation  index,  the  normalized  intensity  variance  of  the  optical  wave,  is  small,  and 
diffractional  effects  must  be  included  (10:894),  (14:21 1 1).  Thus,  whenever  the  scintillation  index. 


,2  _  -  {If 

(/>2 


(2.3) 


where  (•)  denotes  ensemble  averaging,  is  less  than  0.3,  it  is  generally  agreed  that  the  method  of 
smooth  perturbations  is  applicable.  Appendix  B  contains  a  more  complete  discussion  of  the 
method  for  the  interested  reader.  The  accuracy  of  the  smooth  perturbation  method  has  been  called 
into  question  by  some  authors  (2:516).  However,  an  asymptotic  solution  approach,  a  variation  of 
the  method  of  smooth  perturbations,  used  to  solve  Eq.  2.2  in  the  case  of  strong  intensity 
fluctuations  has  become  available  recently  (15:2111).  Additionally,  other  methods  including 
fiiiKtional  integration  methods  have  shown  the  Rytov  method  to  be  quite  accurate  under  certain 
constraints.  These  constraints  are  complex  and  depend  upon  both  the  size  of  the  root-mean-square 
phase  fluctuations  and  the  ratio  of  the  longitudinal  scale  size  of  the  inhomogeneities  to  the  extent 
of  the  Fresnel  zone  of  the  optical  wave  (10:895,896).  The  solution  of  Eq.  2.2  by  the  method  of 
smooth  perturbations  is  still  a  stochastic  quantity  as  the  solution  depends  upon  the  statistical 
characteristics  of  the  random  medium.  It  is  from  this  solution  that  the  statistical  quantities  of 
interest  are  calculated.  It  is  possible  to  bypass  the  solution  of  Eq.  2.2  entirely  and  calculate  the 
statistical  quantities  directly.  This  method  will  be  explored  in  the  following  section. 
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Method  of  Markov  Approximations 

The  method  of  Markov  approximations  facilitates  the  determination  of  closed  partial 
differential  equations  for  the  statistical  quantities  of  the  optical  wave  held.  The  derivation  of 
these  statistical  moment  equations  eliminates  the  necessity  of  explicitly  solving  Eq.  2.2.  As 
shown  in  Appendix  C,  the  statistical  equations  are  obtained  by  approximating  the  optical  wave 
propagation  in  the  inhomogeneous  random  medium  as  a  Markov  process  (29:213),  (25:79).  The 
method  of  Maricov  approximation  is  according  to  (7:1234)  and  (8:171)  a  more  accurate  method  to 
solve  Eq.  2.2  than  the  method  of  smooth  perturbations.  In  order  for  the  method  of  Markov 
approximations  to  be  applicable  to  a  propagation  problem,  it  is  necessary  that, 

1.  The  characteristics  of  the  wave  held  are  dependent  only  upon  the  section  of  the  medium 
through  which  the  wave  has  already  propagated. 

2.  The  longitudinal  correlation  radius  of  7  is  small  compared  to  the  characteristic  length  of 
the  changes  in  ti(p,  z)  (the  complex  amplitude)  in  the  z  direction  (25:79). 

The  hrst  condition  is  actually  impossible  for  the  general  situation  of  optical  wave 
propagation  in  a  random  medium  such  as  the  atmosphere.  Atmospherically  propagating  waves 
will  be  backscattered  to  a  small  degree.  These  backscattered  waves  will  interact  later  with 
forward  propagating  waves  and  therefore,  the  forward  propagating  wave  will  be  dependent  upon 
the  characteristics  of  the  medium  through  which  it  has  not  yet  encountered.  Therefore,  a  direct 
implication  of  the  first  requirement  necessitates  the  consideration  of  solutions  for  Eq.  2.2  as 
opposed  to  Eq.  2.1  (28:5.1 1). 

The  second  requirement  simply  means  that  the  correlation  function 
B({z,  z')  =  (e  (p,  z)  c  (p,  z'))  is  non-zero  over  an  interval  which  is  small  compared  to  the 
characteristic  length  of  changes  in  the  complex  amplitude,  u(p,  z),  of  the  optical  wave  in  the 
mean  direction  of  propagation.  This  in  turn  requires  that  "longitudinal  correlation  radius  of  e  be 
the  smallest  dimension  in  the  problem"  (29:215).  When  this  holds, 

Bt{z,  z')  =  (e  (p,  z)  c  (p,  z'))  can  be  approximated  by  a  scaled  Dirac  delta  distribution,  in 
which  case  7  is  said  to  be  delta  correlated.  Clearly  one  advantage,  which  is  often  exploited  in  the 
literature,  is  the  ability  to  extend  the  integration  range  of  a  function  with  small  support,  such  as  a 
Dirac  delta  fuiKtion,  from  a  finite  to  an  infinite  integration  range.  This  expansion  in  the 
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integration  range  allows  the  approximation  of  certain  correlation  functions  as  Fourier  Transforms 
of  power  spectral  densities  (29:218). 

If  the  previous  two  necessary  conditions  are  approximately  satisfied,  then  the  Maikov 
approximation  method  may  be  applied  to  situations  where  the  scintillation  index  exceeds  unity. 
That  is,  when  the  method  of  smooth  perturbations  fails  to  hold,  namely,  strong  intensity 
fluctuations  (10:894).  The  accuracy  of  the  resulting  equations  of  the  statistical  moments  will  be 
dependent  to  the  degree  that  the  above  two  assumptions  are  satisfied.  Finally,  just  as  the  method 
of  Maikov  approximations  eliminated  the  need  to  explicitly  solve  for  the  solution  to  Eq.  2.2,  the 
application  of  numerical  modeling  methods  and  modem  computer  hardware  has,  in  some 
circumstances,  eliminated  the  need  to  even  consider  the  statistics  of  the  solution. 

Computer  Modeling 

The  application  of  Fast  Fourier  Transform  algorithms,  finite  element  methods  and 
multi-gigaflop  performance  levels  in  modem  computing  hardware  has  allowed  the  computer 
modeling  of  many  random  medium  propagation  problems.  It  is  possible  to  numericaUy  solve 
Eq.  2.1  over  all  locations  in  a  given  volume  simultaneously  using  modem  hardware  and  software. 
However,  even  with  today’s  most  powerful  computing  systems,  it  is  Eq.  2.2  rather  than  Eq.  2.1 
which  is  modeled.  In  one  method,  the  statistical  equations  for  the  optical  wave  are  obtained  in 
integral  form  and  the  method  of  Monte  Carlo  integration  is  applied  to  the  equations  (2:517).  A 
different  method  obtains  solutions  which  are  allowed  to  propagate.  The  solutions  to  Eq.  2.2  can 
be  obtained  relatively  easily  for  arbitrarily  thin  volumes  transverse  to  the  mean  propagation 
direction.  These  solutions  are  then  propagated  to  adjacent  thin  parallel  volumes  sequentially  just 
as  the  actual  optical  wave  would  propagate  through  the  medium  using  known  laws  for  optical 
propagation  (15:21 1 1).  However,  computer  simulation,  though  a  powerful  tool,  has 
disadvantages  over  the  previously  mentioned  analytical  methods.  Simulations,  as  accurate  as 
they  can  be,  often  conceal  fundamentally  important  mechanisms  in  the  propagation  of  optical 
waves.  Additionally,  aliasing,  the  contamination  of  lower  frequency  components  by  higher 
frequency  components,  is  often  an  ignored  problem  when  modeling  continuous  phenomena  using 
discrete  data  on  computers.  Presently  there  is  great  interest  in  one  especially  powerful  method 
from  the  field  of  quantum  mechanics  which  can  be  used  to  model  and  solve  Eq.  2.2  in  a  form 
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suitable  for  numerical  or  asymptotic  analysis  methods  (6:38-1).  This  method,  known  as 
functional  integration,  is  investigated  in  the  next  section. 

Functional  Integration 

Richard  Feynman,  in  his  doctoral  dissertation  in  1942,  proposed  a  novel  solution  to  the 
Schrddinger  wave  equation,  as  shown  in  Appendix  D.  His  solution  allowed  the  number  of 
integrations  in  a  finite-fold  integral  to  diverge  to  infinity  (16:48),  (12:298),  (29:227),  (23:4). 
Feynman’s  technique  spurred  the  investigations  of  others  into  the  area  of  research  which  later 
came  to  be  known  as  "integration  in  a  function  space"  (16:48).  Feynman’s  original  integral 
solution  of  the  Schrodinger  wave  equation  came  to  be  known  as  a  Feynman  integral,  although  it 
was  determined  that  this  integral  was  simply  a  special  case  of  a  Wiener  integral  (5:127).  The 
study  of  Wiener  integrals,  also  known  today  as  path  integrals,  began  with  the  work  of  N.  Wiener 
in  the  1920’s  (16:48).  Because  of  the  functional  similarity  of  the  parabolic  equation  to  the 
Schrddinger  wave  equation,  the  functional  integration  approach  can  be  used  to  solve  the  parabolic 
equation  Eq.  2.2.  Moreover,  it  has  been  shown  that  the  functional  integration  approach  can  be 
applied  to  the  wave  equation  Eq.  2.1  itself  (6:38-3).  It  is  the  method  of  functional  integration 
applied  to  the  parabolic  equation  with  which  this  thesis  will  be  primarily  concerned. 

A  Justification  from  Fourier  Optics  This  section  will  give  an  intuitive,  nonrigorous 
justification  of  the  functional  integral  technique  to  solve  the  parabolic  equation,  Eq.  2.4  in  a 
homogeneous  medium  (7  {p,  z)  =  0).  The  solution  for  the  more  general  case  of  an 
inhomogeneous  medium  will  be  found  in  chapter  3.  The  solution  to  Eq.  2.4 

{2jk^  +  S7l)u{p,z)  =  0  (2.4) 

is  given  by 

where  u{p,  z)  is  the  complex  amplitude  of  the  optical  wave  at  the  z  =  z  plane,  u{p,  0)  is  the 
complex  amplitude  of  the  incident  optical  wave  at  the  z  =  0  plane,  p  =  p(x,  y)  is  a 
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two-dimensional  position  vector  in  the  r  =  ^  plane,  p'  =  is  a  two-dimensional  position 

vector  in  the  z  =  0  plane,  j  =  v/^.  A:  =  x  ^  optical  carrier  wavelength. 

Eq.2.5  is  a  solution  which  is  valid  for  all  positive  values  of  z.  If  one  considers  the 
geometry  in  Fig.  2.1,  then  it  is  easy  to  see  that  if  the  complex  amplitude  of  a  propagating  optical 
wave  is  known  at  the  z  =  0  plane,  then  it  can  be  obtained  at  the  z  =  z  plane  using  Eq.  2.5. 
Furthermore,  Eq.  2.5  must  hold  for  all  values  of  z'  where  0  <  z'  <  z. 


Figure  2.1.  Single  interval  propagation 

Because  Eq.  2.5  holds  for  all  values  of  z'  such  that  0  <  z'  <  z,  it  holds  for  z'  =  ^.  This  simply 
means  that  in  order  to  be  able  to  observe  an  optical  wave  at  the  z  =  z  observation  plane,  it  was 
necessary  to  be  able  to  observe  it  at  some  intermediate  location  along  the  path  of  propagation.  So, 
the  optical  wave  will  propagate  from  the  image  plane  to  the  observation  plane  at  z'  =  |.  Then, 
the  optical  wave  will  propagate  from  this  object  plane  to  the  observation  plane  at  z'  =  z.  Thus, 
referring  to  Fig.  2.2  equation  Eq.  2.5  can  be  applied  to  the  separate  domains:  0  <  z'  <  |  and 
2  <  z'  <  z. 

Therefore,  Eq.  2.5  can  be  applied  to  each  of  the  half-length  sections  and  a  four  fold  integral  is 
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Figure  2.2.  Double  interval  propagation 


obtained 

jj  J J  exp  j  „',^,0)dW^P2 

(2.6) 

where  the  region  of  integration  is  x  R^.  It  is  possible  to  see  the  emergence  of  a  path.  The 
complex  amplitude  at  the  point  (x,  y)  is  determined  by  the  concatenated  rays  or  path 
{pi  —  p),{p2  —  Pi)-  It  is  possible  to  continue  to  divide  the  original  interval  into  three  sections, 
four  sections  and  so  on.  For  the  case  of  an  arbitrary  number  of  sections,  the  geometry  is  shown  in 
Fig.  2.3. 

As  the  number  of  intervals  becomes  large,  the  length  of  each  interval  tends  toward  zero. 
Thus,  the  finite-fold  integral  becomes  a  functional  integral  or  path  integral  in  the  limit  as  the  size 
of  each  interval  vanishes.  The  path  in  this  case  being  any  possible  ray  trajectory  />„,/?„_] , ...,  p\ 
from  the  furthest  observation  plane  backward  to  the  original  object  plane.  The  intermediate 
planes  each  alternately  act  as  object  and  observation  planes  and  in  the  limit,  the  set  of  polygonal 
paths  so  generated  will  almost  span  the  set  of  all  possible  continuous  paths  from  the  original 
object  plane  to  the  final  observation  plane  uniformly.  Formally,  the  solution  to  the  parabolic 
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Figure  2.3.  Multiple  interval  propagation 


equation  is  given  t;'  Eq.  2.7 

"■*(2.7) 

where  the  integration  region  is  lim„_ooR^".  Of  course  for  free  space,  its  easy  to  show  that  the 
preceding  path  integral  is  fully  equivalent  to  the  Green’s  function  from  the  Huygens-Fresnel 
integral.  As  it  was  stated  earlier,  the  previous  path  integral  derivation  showed,  in  a  non-rigorous 
way,  the  applicability  of  the  functional  integral  approach.  It  could  be  argued  that  the 
Huygens-Fresnel  formula  would  give  the  correct  results  over  a  single  propagation  interval 
without  the  need  to  resort  to  the  recursive  subdivisions  that  the  functional  integration  method 
encompassed.  This  is  true;  however,  the  previous  example  was  for  a  deterministic,  homogeneous 
medium.  In  the  case  of  a  randomly  inhomogeneous  medium,  the  functional  integral  approach  to 
solving  the  parabolic  equation,  parallels  t<i<:  heuristic  method  of  considering  thin  slabs  of  the 
medium  to  be  approximated  as  phase  screens.  This  idea  will  be  expanded  further  in  the  following 
chapter. 

Up  until  now,  the  term  complex  amplitude  has  been  used  extensively  to  describe  the 


strength  of  the  propagating  optical  wave.  However,  neither  intensity  detectors,  photographic  film 
nor  most  optically  dependent  biological  organisms  can  see  complex  amplitudes.  All  three  ’’see” 
intensities.  Therefore,  it  is  necessary  to  complicate  matters  by  considering  not  only  the  complex 
amplitude  of  a  propagating  wave  and  all  of  its  statistical  characteristics,  but  also  the  product  of 
the  complex  amplitude  and  its  complex  conjugate  and  the  associated  statistical  characteristics  of 
the  intensity.  Because  the  statistical  characteristics  of  the  randomly  inhomogeneous  medium  will 
further  complicate  matters,  it  is  therefore  necessary  to  consider  the  strength  of  the  fluctuations  in 
the  randomly  inhomogeneous  medium. 

Weak,  Strong  and  Saturated  Regimes  It  has  been  shown  that  for  a  Gaussian  random 
phase  screen  has  the  characteristics  as  shown  in  Fig.  2.4  for  an  incident  unity  amplitude  plane 
wave  (15:21 1 1),  (25:74).  From  this  figure  there  are  three  regions  which  are  of  interest.  The  first 
is  the  region  in  which  the  value  of  0^  increases  from  zero  to  one,  the  weak  fluctuation  regime. 
The  second  region,  if  it  exists,  is  the  region  in  which  0^  is  greater  than  one,  the  strong  fluctuation 
regime.  The  final  region,  if  it  exists  is  where  0^  is  asymptotically  decreasing  to  one,  the 
saturation  regime. 

The  functional  integration  approach  to  solving  the  parabolic  equation  in  a  randomly 
inhomogeneous  medium  is  an  especially  powerful  one  as  the  method  can  be  applied  to  all  three 
regions  of  intensity  fluctuations,  the  weak,  the  strong  and  the  saturated.  This  is  an  especially 
useful  characteristic  as  it  is  possible  to  compare  the  results  from  the  other  mathematical  methods 
to  the  results  from  the  Feynman  integration  method. 

Regions  of  Applicability 

This  section  will  simply  provide  a  tabular  listing  of  the  applicability  constraints  of  the 
various  methods  previously  listed. 
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Scintillation 

Index 


Strong  medium  fluctuations 


Weak  medium  fluctuations 


Distance  from  source 


Weak 

fluctuation 

regime 


Strong 

fluctuation 

regime 


Saturated 

regime 


Figure  2.4.  Scintillation  Index 
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Method 

Regimes 
solution  is  valid 

Advantages 

Disadvantages 

Geometrical  Optics 

Non- 

diffraction 

Very  simple 
method 

Diffractional  effects 
ignored 

Method  of  Smooth 
perturbations 

Weak 

Relatively  simple 
method 

Limited  to  weak 
fluctuations 

Markov  Approximation 

Weak  and 

Strong 

Eliminates  the 
need  to  solve  the 
parabolic  equation 

Unrealistic  assumptions 
possible 

Computer  Modeling 

All 

Eliminates  the  need 
to  calculate  statistics 

Underlying  mechanisms 
are  obscured 

Possible  aliasing 

Feynman  Integration 

All 

Accurate 
and  intuitive 

Infinite-fold  integrals 
difficult  to  analyze 

Table  2.1.  Characteristics  between  various  solution  methods. 


Summary 

This  chapter  has  provided  a  current  literature  review.  This  review  included  several  of  the 
techniques  used  to  solve  the  problem  of  wave  propagation  in  a  randomly  inhomogeneous  medium. 
The  next  chapter  will  explore  the  characteristics  of  path  integration  as  applied  the  this  problem. 
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III.  Feynman  Integration  applied  to  Optical  Wave  Propagation  in  a  Randomly 

Inhomogeneous  Medium 


Introduction 

This  chapter  follows  an  unpublished  paper  by  Dr.  J.  Gozani,  Dr.  V.  Tatarskii  and  Dr.  V. 
Zavarontnyi  (19:1).  Therefore,  in  order  to  help  the  reader,  who  might  become  hopelessly  lost  in 
the  succeeding  sections,  a  brief  overview  of  the  theoretical  development  will  be  provided.  This 
chapter  will  present  a  development  of  path  integration  applied  to  wave  propagation  in  a  randomly 
inhomogeneous  medium  so  that  the  scintillation  index,  the  normalized  intensity  variance,  for  a 
spherical  wave  may  be  calculated.  This  will  be  accomplished  in  three  parts. 

First,  the  existence  of  the  path  integral  solution  will  be  briefly  mentioned  by 
demonstrating  the  similarity  between  the  Schrodinger  wave  equation,  which  has  a  path  integral 
solution,  and  the  parabolic  equation,  with  which  this  thesis  is  concerned. 

Second,  in  order  to  understand  the  development  of  the  Green’s  function  path  integral  for 
the  fourth  statistical  moment  of  the  held,  a  two-fold  path  integral  solution  to  the  two-fold 
stochastic  parabolic  wave  equation  will  be  developed.  This  development  will  be  accomplished  by 
applying  the  Huygens-Fresnel  integral  solution  for  wave  propagation  in  a  deterministic 
homogeneous  medium  to  thin  sections  of  the  randomly  inhomogeneous  propagation  medium.  It 
will  be  shown  that  as  the  number  of  sections  over  which  the  Huygens-Fresnel  integral  is  applied 
increases,  the  resulting  solution  will  become  more  accurate.  In  the  limit  as  the  number  of  sections 
becomes  unbounded,  the  resulting  solution  will  become  a  path  integral  with  stochastic 
components.  A  similar  development  will  yield  a  Green’s  function  path  integral  with  stochastic 
components.  It  will  be  possible  to  separate  the  resulting  solution  into  two  factors.  The  first  factor 
will  represent  the  effects  of  a  completely  deterministic  homogeneous  propagation  medium  and 
the  second  factor  will  represent  the  effects  of  the  stochastic  deviations  of  the  propagation  medium 
about  its  mean  characteristics.  Then,  each  path  in  the  path  integral  will  be  approximated  by  a 
truncated  Fourier-sine  series  using  a  Taylor-functional  series  expansion. 

Third,  using  the  preceding  development  of  the  two-fold  path  integral  for  the  field  as  a 
guide,  the  fourth  statistical  moment  of  the  field  will  be  calculated.  The  stochastic  components  of 
the  resulting  eight-fold  Green’s  function  path  integral  will  be  averaged  over  the  ensemble  of  all 
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paths  from  the  source  point  to  the  observation  point.  After  the  averaging  is  performed,  the 
resulting  equation  will  be  manipulated  in  a  manner  similar  to  the  two-fold  path  integral  up  to  and 
including  the  approximation  of  every  trajectory  by  a  truncated  Fourier-sine  series.  Then  the 
truncated  Fourier-sine  series  will  next  be  approximated  by  the  first  term  in  the  series.  This 
approximation  will  allow  the  eight-fold  path  integral  to  be  approximated  by  an  ordinary 
three-fold  Riemann  integral  which  has  appeared  in  a  similar  form  in  the  current  literature. 

Lastly,  as  it  is  the  scintillation  index,  namely  the  normalized  intensity  variance,  which  is 
ultimately  desired,  it  behooves  the  reader  to  keep  this  ultimate  goal  in  mind  during  the  succeeding 
derivations. 

Parabolic  Equation  and  Schrodinger  Wave  Equation  Similarity 

As  presented  in  Appendix  D,  the  solution  of  the  Schrodinger  wave  equation  may  be  cast 
into  a  ’’sequential  Wiener  integral  with  an  imaginary  variance  term,  namely  a  Feynman  integral 
(16:55).”  The  one-dimensional  Schrodinger  wave  equation  has  the  form, 

where  t)  is  a  wave  function,  V(i,  t)  is  a  potential  energy  function,  and  m,  j  andft  are 
constants  (see  Appendix  D).  Following  a  derivation  similar  to  that  used  to  derive  the 
one-dimensional  Schrodinger  wave  equation  solution,  it  is  possible  to  show  that  there  is  a 
corresponding  two-dimensional  Feynman  integral  solution  to  the  two-dimensional  Schrodinger 
wave  equation.  The  two-dimensional  Schrddinger  wave  equation  has  the  form, 

where  V’(p,  t)  is  a  two-dimensional  wave  function,  V(p,  t)  is  a  two-dimensional  potential  energy 
function,  ^  and  p  =  {x,  y).  Recall  the  form  of  the  parabolic  equation, 

+  ^\  +  k^c  (p,  z)j  «(p,  z)  =  0  (3.3) 
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where  e  (/?,  z)  is  the  fluctuating  part  of  the  medium’s  dielectric  permeability  and  u{p,  z)  is  the 
complex  amplitude  of  the  propagating  wave.  Dividing  both  sides  by  2jk  results  in, 

It  is  easy  to  see  the  similarity  between  Eq.  3.2  and  Eq.  3.4.  Because  of  this  similarity,  it  is  possible 
to  cast  the  solution  of  Eq.  3.4  into  a  Feynman  integral  (6:38-3).  Interestingly  enough,  it  is  possible 
to  obtain  the  solution  to  the  stochastic  Helmholtz  equation,  without  the  parabolic  approximation 
constraints,  as  a  multi-dimensional  Feynman  integral  (6:38-3);  however,  this  thesis  will  not 
investigate  that  topic  other  than  to  mention  it  so  that  the  interested  reader  may  pursue  it  further. 

Following  a  derivation  similar  to  the  one  in  the  chapter  2  subsection  ”A  justiflcation  from 
Fourier  optics”,  it  is  possible  to  solve  the  stochastic  parabolic  equation  using  Feynman 
integration  directly.  Therefore,  it  is  possible  to  obtain  a  solution  to  the  problem  of  wave 
propagation  in  a  randomly  inhomogeneous  medium  by  either  noting  the  similarity  between  the 
Schrodinger  wave  equation  and  the  parabolic  equation  or  by  a  direct  derivation  using  the  notion 
of  an  infinite  number  of  phase  screen  equivalent  propagation  regions.  The  heuristic  derivation  of 
such  a  direct  solution  follows  in  the  next  section. 

Derivation  of  the  Feynman  Integral  Solution  of  the  Stochastic  Parabolic  Equation 

This  section  will  present  a  non-rigorous  derivation  of  the  solution  to  the  stochastic  (e?^  0) 
parabolic  equation  using  Feynman  integration.  This  derivation  will  establish  the  earlier  claim  of 
an  intuitively  pleasing  solution.  The  interested  reader  can  find  a  more  mathematically  rigorous 
derivation  of  the  relationship  between  the  stochastic  parabolic  equation  and  its  Feynman  integral 
solution  in  (29:227,228),  or  (7:1225-1227).  Moreover,  there  is  a  derivation  and  so^’ition  of  a 
related  acoustical  stochastic  parabolic  equation  available  in  the  literature  (23:6-12).  Lastly,  there 
exists  an  extensive  reference  table  which  relates  quantum  mechanical  functionals  to  acoustical 
functionals  which  appear  in  Feynman  integral  solutions  to  the  Schrodinger  wave  equation  and 
acoustical  parabolic  equation  respectively  (23:19).  As  it  was  stated  in  chapter  2,  the  functional 
similarities  of  the  equations  in  acoustics,  quantum  mechanics  and  optical  wave  propagation  has 
lead  to  similar  general  solutions.  Therefore,  it  is  important  to  consider  the  techniques  and  results 
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from  other  research  fields  for  novel  methods  to  obtain  solutions  to  wave  propagation  in  a  random 
medium.  This  is  why  the  work  in  other  research  fields  has  been  brought  to  the  reader’s  attention. 

Beginning  with  Eq.  3.5  which  is  derived  in  Appendix  E,  the  complex  amplitude,  u(p,  z), 
of  an  initially  unit  amplitude  plane  wave  propagating  in  the  z  direction,  for  a  relatively  short 
distance  through  the  randomly  inhomogeneous  medium,  is  given  approximately  by, 

~  ^ 

where  p"  is  a  vector  representing  the  observation  or  measurement  location  in  the  observation 
plane,  //  is  a  vector  representing  a  point-source  location  in  the  image  plane,  z  is  the  propagation 
distance.  Implicit  in  the  previous  equation  was  the  assumption  that  7  (p,  z)  was  changing  slowly 
across  a  small  propagation  distance  z.  As  the  distance  between  the  image  plane  and  observation 
plane  decreases,  the  above  approximation  will  become  more  accurate.  A  closer  examination  of 
the  above  equation  will  yield  an  interesting  interpretation.  The  effects  of  the  inhomogeneous 
medium  on  the  propagating  unit  amplitude  plane  wave  may  be  duplicated  by  an  equivalent 
idealized  random  phase  screen  placed  before  a  homogeneous  medium,  so  long  as  the  propagation 
distance  is  sufficiently  small  (29:228).  The  actual  dimensions  of  ’’sufficiently  small”  will  be 
explained  later.  An  idealized  phase  screen  is  a  mathematical  construction  which  has  the  physical 
properties  of  vanishingly  small  (differential)  thickness,  and  a  transmittance  characteristic  which 
changes  only  the  phase,  and  not  the  amplitude,  of  an  incident  optical  wave.  Thus,  all  the  effects 
of  the  randomly  inhomogeneous  medium  may  be  modeled  as  having  been  induced  by  an 
equivalent  random  phase  screen.  Therefore,  by  recognizing  this  interpretation,  Eq.  3.5  may  be 
rewritten  with  the  phase  screen  interpretation  as, 

' - - - ' 

Phase  screen  or 
random  phase  plane  wave 

We  can  rewrite  this  equation  in  another  form  by  first  noting  that  for  a  thin  slab  of  the  propagation 
medium  we  have 

c  (p',z)  ~c  (p',0)  (3.7) 
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as  we’ve  assumed  that  e  varies  slowly  in  the  z  direction.  This  being  the  case,  we  can  write 


J I  «(/>'.  0)  exp  y  -  p'?^  (3-8) 

where  exp  c  {p',  z)j  =  u(/y,  0).  Furthermore,  if  a  subdivision  process,  similar  to  the  one 
applied  in  chapter  2  to  the  problem  of  propagation  in  a  homogeneous  medium,  is  applied  to  this 
randomly  inhomogeneous  medium  (Eq.  3.5),  we  will  obtain  a  ’’sequential  >\^ener  integral  with 
imaginary  variance  parameter”  (5:127).  This  is  another  way  of  saying  a  Feynman  integral  will 
result  as  a  solution.  This  solution  is  written  formally  as. 


u{p",z)  =  lim  u{pn,z)  = 


lim 

n— >00 


n— »oo 
nk 


\j2Trz 


exp 


drpm 


(3.9) 


where  p"  has  been  defined  as  lim„_oo  Pn  and  a  more  general  boundary  condition  u(po»  0)  has 
been  included. 


By  interpreting  the  previous  infinite-fold  integral  as  it  was  done  for  Eq.  3.5,  then  it 
becomes  clear  that  the  effects  of  a  randomly  inhomogeneous  medium  on  a  propagating  plane 
wave  will  be  indistin!>uishable  from  the  effects  of  an  equivalent  (though  infinite)  set  of  ideal 
random  phase  screens.  Additionally,  it  was  assumed  that  e  could  be  approximated  as  being 
nearly  constant  across  the  propagation  distance  and  as  the  propagation  subdivision  width 
decreases  and  eventually  becomes  a  differential  in  thickness,  the  approximate  solution,  which 
corresponds  to  an  infinite  number  of  small  propagation  regions,  becomes  an  exact  solution. 
Therefore,  unlike  the  case  of  an  homogeneous  medium,  it  is  necessary  to  let  the  number  of 
subdivisions  and  therefore  the  number  of  equivalent  phase  screens  become  unbounded  to  obtain 
an  exact  solution.  Thus,  it  is  seen  that  the  complex  amplitude  of  an  optical  wave  propagating  in 
an  inhomogeneous  medium  is  represented  as  a  two-dimensional  Feynman  integral.  This  integral 
is,  as  we  have  seen,  is  the  limiting  case  of  replacing  the  random  inhomogeneous  medium  by  one 
containing  a  set  of  uniformly  spaced  random  phase  screens  and  allowing  the  number  of  phase 
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screens  to  become  unbounded. 


The  Feynman  Integral  Solution  in  the  Umit  We  can  make  an  additional  interpretation 
concerning  Eq.  3.9.  This  interpretation  will  result  in  a  functional  integral  form  which  has 
appeared  in  the  literature  (5:135),  (22:1264),  (1:37),  (21:8),  (29:228),  (23:11).  By  considering 
Eq.  3.9  it  follows  that  the  functional  integrand. 


exp 


t  «p  (f  (p™.  v)) 


may  be  rearranged  to  yield, 


exp  (  f  t 

^  m=l 


\Pm  Pm— li 


~  z  z 

e  (pm,m-)- 

n  n , 


(3.11) 


Now,  consider  the  following  interpretation  of  the  above  expression.  Eq.  3.11  is  a  functional  of  the 
n  +  1  variables  pq,p\-’-  Pn-  Equivalently,  it  is  possible  to  consider  Eq.  3.1 1  as  a  functional  of  a 
differentiable  parametric  function  pij)  where  r  is  a  continuous  parameter  evaluated  at  discrete 
locations, 

{POi  P\i  P2i  •••)  Pn)  —  {p('?■)lT=0,r=i,r=2i•••,T=(n-l)i,T=^}  (3-12) 

That  is,  p(r)  is  a  continuous  variable  of  a  continuous  parameter  r  sampled  at  the  evenly  spaced 
points  r  =  0,  T  =  ^,  T  =  2^  •  •  •.  In  the  limit  as  n  becomes  unbounded  in  Eq.  3.1 1,  the  finite, 
discrete  sampling  of  the  functional  p(t)  will  become  effectively  continuous.  The  astute  reader 
will  immediately  disagree  as  to  the  sampling  being  continuous  as,  for  example,  there  will  be  no 
sampling  at  any  of  the  irrational  numbers  due  to  the  limiting  process  described  above.  That  is 
true;  however,  because  any  irrational  number  may  be  approximated  as  closely  as  desired  by  a 
convergent  sequence  of  rational  numbers,  p(t)  will  be  sampled  effectively  everywhere. 

Now  consider  a  simple  substitution  of  6t  for  ^  into  Eq.  3.1 1.  This  will  yield  a  more 
apparent  result  from  elementary  calculus, 

®*P  H  ^  ^  (/Jm,m5T)«T^  (3.13) 

\  fTl=l  /  \  TTI=0  / 
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To  begin,  consider  the  first  exponential  factor  in  the  previous  equation.  It  follows  that 


f  |Pm  Pm-l|\  ^ 

/»(T)|r=JM  -  p(T)|^_(m-l). 

2 

dpir)  ^ 

\  St  J  ~ 

—  t| 

"  n 

dr 

(3.14) 


Therefore,  in  the  limit  as  n  becomes  unbounded  it  follows  that. 


exp 


\Pm  Pm— 1 1 

£ 

n 


(3.15) 


Moreover,  the  second  exponential  factor  in  Eq.  3.13,  may  be  manipulated  in  a  similar  fashion. 
Therefore,  in  the  limit  as  n  becomes  unbounded  it  follows  that, 

exp  ^  ^  Jo  ^  (pW, 'r)dT^  (3.16) 


Therefore,  the  entire  functional  integrand  may  be  written  as, 

where  t  is  the  limiting  continuous  variable  p.  Therefore,  it  is  possible  to  write  Eq.  3.9  formally  as, 

«(P>2)  =  Jo  ^ 


where  I>^[p(t)]  =  lim„_oo  nm=i  This  is  one  of  the  formal  forms  for  the 

Feynman  integral  solution  that  is  often  found  in  the  literature.  It  should  be  noted  that  when  we 
have  a  deterministic  medium,  7=0  and  the  preceding  integral  will  reduce  to  Eq.  2.7.  The  reader 
should  be  aware  that  similar  to  the  Riemann  integral  where  the  integration  variable  is  simply  a 
dummy  variable,  p(t)  is  similarly  a  dummy  integration  variable.  Therefore,  Eq.  3.18  is  fully 
equivalent  to 


«(,,z)  =  //i)^Wexp  f  l^f  exp  (f  I'  ?  (x,r)dr)  (3.19) 
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where  x  =  x(t)  and  t  is  the  limiting  continuous  variable. 


The  preceding  formal  equation,  it  must  be  understood  is  only  that,  formal.  It  has  been 
found  that  the  Feynman  integral  does  not  possess  a  countably  additive  measure  (26:8).  Therefore, 
it  would  be  more  accurate  if  the  preceding  equation  was  rewritten  as 


u(p,  z) 


dx 

dr 


dr^  exp  (y 


(3.20) 


where  "paths”  means  add  all  the  contributions  from  all  paths. 


Green ’s  Function  Path  Integral  Following  the  derivation  of  the  preceding  equation,  it  is 
possible  to  show  that  the  Green’s  function  for  propagation  in  a  randomly  inhomogeneous  medium 
is  given  by 

G(/,2;/)',0)  =  e  (a:(T),T)dT| 

°  (3.21) 

where  each  trajectory  [x]  begins  at  p'  when  t  =  0  and  ends  at  p"  when  t  =  z,  the  overdot 
indicates  differentiation  with  respect  to  r  and  each  trajectory  was  the  limiting  case  of  a  polygonal 
curve.  The  total  complex  field  amplitude  at  an  observation  point  due  to  several  point  sources  or  to 
a  continuum  of  point  sources  in  a  source  plane  some  distance  away  is  then  given  by 

u{p",z)  =  11  G{p'\  z;  p',  0)u(p',  0)dV'  (3.22) 


which  is  simply  a  superposition  integral  with  a  Green’s  function  integrand.  Let  the  trajectory 
x(t)  be  decomposed  into  two  components,  x(t)  and  x  (t).  The  first  component,  x(t)  is  the 
straight  line  connecting  a  given  source  point,  pf  with  the  observation  point  p".  The  equation  of 
this  line  is  simply 


x(r)  = 


(3.23) 


where  t  £  [0,  z].  The  second  component,  x  (r)  is  simply  the  difference  between  x(t)  and  x(t). 
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that  is,  the  fluctuating  component.  Therefore,  a  given  trajectory  x(t)  may  be  written  as 


x(t)  =  i(t)+  X  (t) 


(3.24) 


It  is  important  to  notice  that  the  fluctuating  component  x  (r)  vanishes  at  the  boundaries  r  =  0 
and  T  =  z.  It  is  important  to  notice  that  the  strait  line  component  corresponds  to  a  completely 
deterministic  propagation  medium  and  the  fluctuating  component  corresponds  to  the  stochastic 
deviations  of  the  propagation  medium  about  its  mean  characteristics.  Applying  this  trajectory 
decomposition  to  every  trajectory  in  the  Green’s  function  in  Eq.  3.21,  we  obtain 


G{p",z-,p',o)= 

J  Jx(Q)=Q 

+2J(t)- 

exp{y  €  (i(t)-I-x  (T),T)dr| 

(‘^)H 

(3.25) 


where  we  have  used 


ii(T)p  =  i(T)-x(T) 

=  ^^(r)  +  x(t)^  .  ^j(T)  +  x(t) j 

=  |^(r)p  +  |x(T)p  +  2^(T)-x(T) 

'W'-pV"" 


(3.26) 


Let’s  consider  the  first  exponential  term  in  the  path  integral  above. 


dr 


(3.27) 


We  can  simplify  this  equation  by  first  noting 

0*:  rf\p"-p'\^^ 


exp 


(3.28) 
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and  second  noting  that 


=  exp(O)  =  1 


(3.29) 


The  integral  in  the  preceding  exponential  vanishes  due  to  the  boundary  conditions  on  x  {t).  We 
can  now  write  the  Green’s  function  path  integral  as 

Gip",  z; p', 0)  =  exp  -  pf )  j D^[x]  (3.30) 

(jk  /"^  |~  P  ~  ~ 

~2  Jo  +  «  (x(t)+i  (r),T)dT 

For  free  space  (e=  0),  the  Green’s  function  is  known  to  be 

Goip",z-,p\0)  =  ^exp  (I^Ip"  -  pf )  (3.31) 


where  the  subscript  0  is  used  to  denote  free  space.  Therefore,  we  can  write  the  Green’s  function 
for  a  general  medium  as 

G(p",z;p',0)  =  Go(p",z;p',0)F(z)  (3.32) 


where 


f  D^[^]  exp 

K  J  Jx(0)=0 


+  7  (i(t)+  I  (T),T)dr 


1 


(3.33) 


We  will  be  interested  in  the  statistical  characteristics  of  F{z)  because  Go{p",  z\ p',  0)  is 
completely  deterministic. 

As  an  aside,  the  astute  reader  will  immediately  realize  that  the  free  space  Greens  function 
will  supply  a  normalizing  condition  on  the  path  integral  of  Eq.  3.30.  Specifically,  when  the 
propagation  medium  is  free  space,  we  have  €  =  0.  Therefore,  it  follows  that  for  free  space 


Go(p",2;p,0)  =  exp  ^^|p" -p'l^j  /.^(o)-o  1^^^^  } 


(3.34) 
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which  further  implies  that 


r  r®(i)=o 

J  Jx(0)=Q 


D^[a:]exp 


1 

jXz 


(3.35) 


Before  we  begin  the  next  section,  it  will  prove  enlightening  to  verify  that  the  Green’s 
function  path  integral  will,  at  least  for  the  case  of  free  space,  reduce  to  the  regular  Green’s 
fimction  which  is  derivable  from  the  Huygens-Fresnel  integral.  Additionally,  it  is  hoped  that  this 
derivation  will  help  convince  the  reader  that  the  paraxial  approximation  upon  which  the 
Huygens-Fresnel  integral  is  based  and  the  paraxial  approximation  upon  which  the  Feynman 
integral  is  based  are  equivalent.  To  begin,  we  know  that  the  Green’s  function  for  the 
Huygens-Fresnel  integral  is  given  by 


Ghf{p",z,p\0)  = 

where  the  HF  subscript  will  be  used  to  denote  the  ’’Huygens-Fresnel”  Green’s  function.  Now, 
we  will  form  the  sequence  of  approximations  to  the  Green’s  function  and  then  show  that  each 
approximation  is  in  fact  equal  to  every  other  approximation.  Specifically,  consider  the  case  where 
there  is  one  subdivision  for  the  trajectories.  In  this  case,  let  /oq  be  a  general  location  vector  in  the 
source  plane  and  be  a  general  location  vector  in  the  measurement  plane.  As  shown  in  chapter 
2,  we  have  for  the  complex  field  amplitude, 

u{p\,z)  =  j I u{po,0)txp  (^^\p\  -  d^Po  (3-37) 

By  using  the  property  of  the  Green’s  function, 

«(/,  z)  =  Jj  Gip^',  z,  p\  0)u(p\  0)dV  (3.38) 

we  see  that  the  resulting  Green’s  function  is  equal  to  Ghf-  Now  consider  a  propagation  region 
with  two  subdivisions.  In  this  case  let  /tq  be  as  before  a  location  vector  in  the  source  plane. 
Additionally,  let  pz  replace  p\  as  a  location  vector  in  the  measurement  plane.  Lasdy,  let  pi  now 
represent  a  general  location  point  which  is  midway  along  the  propagation  medium  between  po 
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and  p2.  In  this  case,  the  complex  field  amplitude  in  the  measurement  plane  which  will  be 
dependent  upon  the  complex  field  amplitude  at  the  intermediate  (pi)  plane  will  be 


«(/>2,«)  = 

^1  ) 


it  follows  that  the  Green’s  function  may  be  written  as 

2 


g(/>2,2,a),o)  =  + 


Now,  consider  the  factor 


which  may  be  rewritten  as 


Pi 


|p2  -piP  +  |pi  -  Pop 


When  the  preceding  factor  is  substituted  into  Eq.  3.40,  we  obtain 

2 


G(«,  z,  A,,  0)  =  (^)  exp  -  pop) 

//'•“'■  V2x/2 


expf^|p.-(^')n^V, 


which  is  a  gaussian  integral.  Such  an  integral  is  easily  evaluated  and  yields 


(3.41) 


(3.42) 


(3.43) 


G(p2,2,Po,0)  =  ^^^exp  ^i^|p2  -  poP^  =  GffF(p2,^',Po,0)  (3.44) 


the  free  space  Green’s  function.  Lastly,  it  is  relatively  simple,  though  extremely  tedious,  to  show 
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that  for  N  subdivisions,  the  Green’s  function 


G(W,«,A),0)  = 

//  ■  ■  ■  //“"  (^  ^ I')  n 


(3.45) 


can  also  be  simplified  so  that 


G{PN,  z,PQ,0)  —  Ghf{pn,  2,  A)>  0)  (3.46) 

Therefore,  because  the  choice  of  N  was  arbitrary,  we  can  take  the  limit  as  N  becomes  unbounded 
and  we  conclude  that 

G(p",z;p',0)  =exp^^|p"-p'|2^ 
therefore, 

\i{r)f  dr  j  ^  (3.48) 

and  so  we  conclude  that  the  Green’s  function  from  the  Huygens-Fresnel  integral  and  the  Green’s 
function  from  the  Feynman  integral  are  equal  for  free  space. 

Orthogonal  Function  Expansion  It  is  known  that  if  a  function  satisfies  certain 
requirements,  then  it  possible  to  expand  it  over  an  infinite  orthogonal  function  set.  The  choice  of 
the  orthogonal  functions  and  therefore,  whether  the  set  is  countable  or  not,  is  usually  chosen  to 
simplify  the  problem  at  hand.  For  our  problem,  we  will  choose  a  Fourier-sine  series  as  the 
orthogonal  function  set.  This  set  will  be  used  to  perform  an  expansion  of  the  fluctuating  part  of 
each  trajectory  in  the  Green’s  function  path  integral. 

Let  the  fluctuating  part  of  the  trajectory  x  (t)  be  expanded  in  a  Fourier-sine  series  as 

^  ("t)  = 

n=l 
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where 


and  Xn  GR^.  The  use  of  the  above  Fourier-sine  series  results  in  the  simplification 

expjYjf  |i(r)|  drj  =  exp|i^Ji„p| 

as  the  cross  terms,  all  being  orthogonal,  will  integrate  to  zero. 


(3.50) 


(3.51) 


Approximation  of  the  Green’s  Function  Path  Integral  In  this  section,  we  will  use  a  finite 
truncation  of  the  Fourier-sine  series  to  obtain  an  approximate  value  of  z;  0).  Consider 
the  path  integral 


f  («)  =  f  exp  "T  /  ^  ^3-^2) 

JJx(o)=o  2  Jo  } 


where 


y 

(t)  =  52a:n'^n(T) 

n=5l 


(3.53) 


We  seek  to  find  a  finite  value  of  N  such  that  Iqo  (t)  can  be  replaced  by  ijv  (t)  in  Eq.  3.52.  This 
is  not  a  simple  problem  as  the  choice  of  N  will  be  dictated  by  several  factors  including  the 
statistics  of  e.  This  problem  not  withstanding,  we  may  make  some  simplifications  to  Eq.  3.52. 

Let  us  begin  with  the  integral  term  in  the  exponent  of  the  path  integrand  in  Eq.  3.52 


i: 


€  ix{T)+  Xqq  (T),T)dT 


(3.54) 


It  is  seen  that  e  is  a  functional  of  the  functions  x  and  xat.  Therefore,  by  following  the  method 
given  in  Appendix  F  for  expanding  an  integral  of  a  functional  in  a  variational  Taylor  series  we  get 
for  the  zeroth  variational  term 


/o' 


7  (l(T)-f  I/v  (T),T)dT 


(3.55) 
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while  tile  first  variational  term  is 


a  €  (l(T)+lAf 
dxN  (t) 


(3.56) 


and  the  second  variational  term  is 


1  r  ~  ~  (a:(T)+  (t),t)  ^ 

2!  Jo  <"> - P-''” 


where 

X-^  (t)  =  Xn<j>n{x)  (3.58) 

n=JV+l 

Higher  order  terms  are  similarly  defined  as  the  interested  reader  will  find  in  Appendix  F.  Similar 
to  the  more  familiar  Taylor  series  in  which  the  value  of  a  function  at  a  fixed  point  is  approximated 
by  the  first  term  of  the  series,  it  is  possible  to  approximate  E(].  3.54  by  the  first  term  of  its 
variational  series  given  by  Eq.  3.55.  Additionally  the  other  exponential  factor  in  the  path 
integrand  of  Eq.  3.52  may  be  simply  truncated  after  N  terms.  Doing  so  results  in 

-  f  L  D^^n]  exp  ( I  51 1®"!^  +  "t  /  ^  (r),  r)dT|  (3.59) 

jJxn{o)=o  2  yo  J 


This  result  will  be  used  in  the  next  section  where  this  thesis  will  be  concerned  with  the  fourth 
statistical  moment  of  the  complex  amplitude  of  a  propagating  optical  wave.  The  ultimate 
accuracy  of  these  approximations  will  be  determined  by  comparing  the  predictions  of  these 
approximations  for  scintillation  index  with  reality.  Before  we  investigate  the  fourth  statistical 
moment  of  the  held,  it  should  be  noted  that  the  previous  equation  is  the  result  of  making  two 
different,  though  related  types  of  approximations.  The  first  approximation  was  made  by  assuming 
every  trajectory  of  the  Feynman  integral  could  be  approximated  by  a  truncated  Fourier-sine 
series.  The  second  approximation  was  made  by  assuming  that  the  difference  between  the  € 
functional  of  the  original  trajectory  and  the  e  functional  of  the  approximate  trajectory  was 
sufficiently  small  so  that  the  former  could  be  approximated  by  the  later.  In  other  words,  the 
approximation  by  the  first  term  in  the  variational  Taylor  series  was  sufficiently  accurate. 
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Fourth  Statistical  Moment 


The  fourth  statistical  moment  of  the  held,  which  gives  the  intensity  variance  of  the 
propagating  optical  wave,  will  be  investigated  in  this  section.  The  other  moments  of  major 
interest,  the  first  and  second,  can  be  obtained  in  a  manner  similar  to  that  which  will  be  used  to 
obtain  the  fourth  moments.  Therefore,  this  thesis  will  consider  only  the  derivation  of  the  fourth 
moment  of  the  held.  The  derivation  of  the  fourth  statistical  moment  of  the  held  will  culminate  in 
the  calculation  of  the  scintillation  index  for  a  point  source.  Thus,  the  general  fourth  statistical 
moment  will  be  obtained  and  then  the  four  independent  points  in  the  general  fourth  moment 
equation  will  be  meiged  into  a  single  point.  This  is  the  ultimate  goal  of  which  the  reader  must  be 
mindful  during  the  reading  of  the  succeeding  sections. 

It  is  reasonable  to  wonder  why  the  general  fourth  moment  equation  is  obtained  and  then 
the  four  independent  points  replaced  with  a  single  point  rather  than  begin  with  a  single  point  in 
the  hrst  place.  After  all,  the  point  source  scintillation  index  or  normalized  intensity  variance,  is 
ultimately  desired.  The  necessity  of  doing  so  for  the  extended  medium  parallels  the  necessity  of 
doing  so  for  a  phase  screen.  The  closed  differential  equation  which  describes  the  fourth  statistical 
moment  for  the  phase  screen  case  is  dependent  upon  the  second  partial  derivatives  of  the  four 
independent  points.  Is  is  only  after  the  appropriate  equations  have  been  solved  that  in  the  limit, 
the  four  points  are  allowed  approach  eaci,  other  and  ultimately  be  replaced  by  a  single  point 
(25:22). 

As  we  have  already  seen,  the  field  at  an  observation  point  may  be  written  in  terms  of  all 
point  source  contributions  in  the  object  plane.  Specifically  this  was  a  superposition  integral  of  the 
Green’s  function  of  the  problem.  Therefore,  the  fourth  moment  of  the  field  may  be  written  as 

r4  =  {Uip'{,z)U*{f/[,z)U{p'i,z)U*{plz))  (3.60) 

where  (■)  denotes  ensemble  averaging.  Assuming  the  mean  direction  of  propagation  is  the  z 
direction,  we  have  for  the  field 


U{p,z)  =  u{p,  z)  exp(jkz) 


(3.61) 
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Using  this  substitution  we  obtain 


TAip'i  pI  pI  pI  z)  =  z)vr{p'i,  z)u{(/i,  z)u-{p'i  z))  (3.62) 

This  may  be  written  in  terms  of  the  appropriate  Green’s  fimctions  and  superposition  integrals 
resulting  in 

^A{PuP'i,p'lp'iz)^  (3.63) 

II II II  !!^G{pU'^p\,Q)G*{f/i,z^fi,,Q)G{pU^^ 

«(Pi,0)u*(/4,0)«(p'3,0)u*(p;,0)(iVi<iV2dV3d^p; 

Let  the  statistical  moment  of  the  Green’s  functions  be  denoted  as 


G<4>  —  G<4>  (Pj  ,  P2^  P3  i(^At^‘*  p\ » P2i  Pv  p'ai  0)  = 

(G(pr,z;p'i,0)G-(p^,z;pi,0)G(p?,2;p^,0)G*(p:(,z;p;,0)) 


(3.64) 


Each  of  the  Green’s  functions  in  the  above  ensemble  average  may  be  written  in  terms  of  the 
appropriate  path  integral  which  results  in  the  fourth-order  path  integral 


_  f  rxM=p"  f  r  ^x,(z)=p"  f  [H(z)=p'l 

G<4>  =11  I  II  11  D^[xi]D^[x2]D^[x-i]D^[x4]  (3.65) 

e*p  |y  (|*i(t)1^  -  |i2(T)P  +  li3(T)p  -  li4(T)P)  dr^ 

^exp  |y  e  (ii(t),  t)-  c  (x2(t),  t)+  e  (i3(t),  t)-  c  (x4(t),  T)(iT|^ 


If  for  the  moment,  let  us  consider  the  stochastic  quantity  E 


^  ^  ^  '  (^3(t),t)-C  (x4(T),T)dT|^  (3. 


66) 


Let  us  assume  that  the  Markov  approximation,  as  described  more  fully  in  Appendix  C  holds  for 
6.  Let  us  also  assume  that  7  is  a  zero  mean  gaussian  random  held.  Then  the  above  ensemble 
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average  may  be  simplified.  Specifically,  we  assume  a  correlation  function  of  the  form 


5e(p  -  p',z-  z’)  =  (1  (p,z)  €  (p',z'))  =  A{p  -  p')6{z  -  z')  (3.67) 


where  is  the  Dirac  delta  distribution.  Therefore,  we  assume  that  e  is  a  gaussian 
delta-correlated  random  field  as  in  (29:214).  The  gaussian  assumption  is  a  good  assumption  as 
the  effects  of  each  small  section  of  the  randomly  inhomogeneous  medium  can  be  modeled  as  an 
independent  identically  distributed  random  variable.  The  assumption  of  delta  correlation  is  a 
reasonable  first  approximation  as  it  implicitly  assumes  no  backscattering  effects.  These 
assumptions  allow  the  simplification  of  the  statistical  quantity  E 

E  =  exp|-|-  (c  (a:i(T),T)-  e  (a:2(T),T)+ e  e 

Expanding  terms,  integrating  and  taking  averages,  results  in 

£;  =  exp^-^  F[xx{T),X2iT),XiiT),x4T),T]dT 


where 


F[ii  (t),  i2('r),  a:3(’‘),  M^),  r]  =  4.4(0)-  (3.70) 

A{xx(t)  -  X2{t))  -I-  A{xx{t)  -  X2(t))  - 
.4(ii(t)  -  X4(t))  -  Aix2{T)  -  i3('’r))  + 

A{x2{t)  -  X4(t))  -  Aix^ir)  -  X4(t)) 


and 

>1(0)  -  A(p)  =  cos(jk  •  p))^,{k,0)d^k  (3.71) 

where  <>c(A;,  0)  is  the  two-dimensional  power  spectral  density  of  the  the  fluctuations  in  e. 
Furthermore,  for  a  pure  Kolmogorov  law 


>l(0)-.4(p)  =  1.46CfHf 


(3.72) 
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under  the  Markov  approximation  (29:219).  Returning  to  the  fourth  order  path  integral  we  have. 


r  r  rx2(z)=i^  t  r 

J  .'ii(0)=p',  J  Jx2(0)=fifj  J  Jxi{0)=p'j  J  /e4(0)=p' 

D^[xi]D^[x2]D^[x^]D^[x4] 

exp  |y  jf  |ii(t)|2  -  |a:2(T)|^  +  |a:3(T)P  -  |x4(T)pdT| 

{  Jo 


(3.73) 


exp ' 


Examination  of  the  previous  equation  will  reveal  that  the  first  exponential  term  which  is 
independent  of  the  fluctuations  of  the  medium  represents  a  pure  diffraction  factor.  This  would  be 
the  only  factor  if  propagation  were  occurring  in  free  space.  The  second  exponential  factor 
represents  combined  effects  of  both  diffraction  and  refraction.  This  term  is  the  scattering  function 
of  the  fourth  moment  along  trajectories  influenced  by  refraction  and  diffraction.  Furthermore,  for 
short  distances,  such  as  the  differential  width  of  a  phase  screen,  the  diffractional  effects  of  the  pure 
diffractional  factor  dominate  the  diffraction-refraction  factor.  For  longer  propagation  distances, 
this  will  in  general  no  longer  be  true;  however,  it  is  not  correct  to  state  that  at  long  propagation 
distances  the  combination  diffraction-refraction  (scattering)  term  will  dominate  the  pure 
diffraction  term,  but  rather  that  the  pure  diffractional  factor  will  not  necessarily  dominate  (18:1). 

The  previous  equation  which  represents  the  statistical  effects  of  a  propagation  medium  on 
four  beams,  may  be  simplified  by  the  following  variable  substitution 


4jR(r)  =  X|(t)  +  X2{t)  -I-  Z3('’')  +  X4{t)  (3.74) 

P(t)  =  h(t)  -  X2(t)  +  X3(t)  -  X4(t) 

2ri(T)  =  Ii(t)  +  X2(t)  -  X3(t)  -  X4(t) 

2r2(T)  =  xi(t)  -  X2(t)  -  X3(t)  +  X4(t) 


These  4  substitutions  result  in 

|xi(t)P-|x2(t)P-|-|x3(t)P-|x4(t)P  =  2.R(t)  • /j(T)  +  2rj  (t)  •  r2(T)  (3.75) 
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and 


F[xiiT),X2{T),x^(T),X4iT),T]  =  F[ri  (t),  r2(T),  p(t),  t]  =  (3.76) 

^(ri(T)  +  ^p{r))  +  i4(r,(T)  -  ^p{t))  + 

Air2iT)  +  ^p(t))  +  A{r2{T)  -  |/)(t))  + 

.4(ri(T) +  r2(T))  -  i4(ri(T)  -  r2(T)) 


It  is  seen  that  the  scattering  factor  will  be  independent  of  R{t)  under  the  preceding  variable 
transformation.  At  the  present  stage  of  development,  the  Green’s  function  is 


r  fR(z)=K'  r  rp(z)=p-  r  /‘riU)=r-,-  e  fri(z)=r-^ 

G<4>  =  //  //  //  I  D^[R]D^[p]D^[r^]D^[r2]  (3.77) 

JJlH0)=R'  J  Jp(Q)=pf  J  Jrtm=r',  J  Jr2(0)=r', 


|■P(^)=P" 

Jp(0)=p’ 


^r,(z)=r'/ 

/r,(0) 


r2(z)=r2 


exp 


••2(0)= 


R{t)  •  p(t)  +  ri(T)  •  r2(T)  -  —F[ri{T),r2{T),p{T),T]dT 


I 


As  it  was  done  previously,  consider  the  trajectories  R{t)  and  p(r)  as  strait  lines 
connecting  source  and  observation  points  with  some  deviation  superimposed  upon  them.  Then 
we  can  make  the  following  substitutions 


R{t) 

«(t) 

p{t) 

p(t) 


=  R{t)+  R  (t) 
= 

=  p(t)+?(t) 


(3.78) 


where  R!'  =  i2(T)l^_^ ,  iZ'  =  /i(T)|^_(,,  and  so  fourth.  Using  these  substitutions,  it  follows 
directly  that 
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exp{jA:jr  ^(t)  •  p(T)dT|  = 

exp  •  ip"  -  p')  J  exp  jjT:  jf  Rir)  •  P(T)dT| 


(3.79) 


as  the  cross  terms  will  integrate  to  zero.  Furthermore,  the  remaining  integral  may  be  integrated 
by  parts  which  results  in 

expjjifcjf  ^(t)  •  P(T)dT|  =  exp  j-jJA:  jT  S(T)-p(T)dT|  (3.80) 


The  Green’s  function  now  takes  the  form 


G<4>  =  exp  |^(iZ"  -  m)  •  ip"  -  p')}  (3.81) 

/■fi{i)=K"  fP(z)=p"  f  fr\(z)=r\'  r  ^r2(z)=r''  ^  , 

h  ~  //  I  .  D^[R]D^[p]D\,]D^[r2] 

Jr(0)=IV  J P{Q)=pf  J  JT\(Q)=r’^  J  7rj(0)=rJ 

exp  j-i*:  R  (t)  ■  Pi'^)  -  ni'^)  •  r2(T)(iTjexp  ^  i^[ri(T),r2(T),p  (t)  +  p(T)]dT| 

A  path  integral  simplifying  formula  referred  to  as  the  ’’delta  functional”  property  of  path  integrals 
exists  which  can  be  used  to  further  simplify  the  preceding  path  integral  (9:22).  This  path  integral 
relation  is  given  by 

j I  D[xiz)]D[yiz)]  exp  jf  1(2)  •  (y(z)  -  fiz))dz^  F[y{z)]  =  i^[5(«)] 


(3.82) 


where  1(2),  y(2),  /{z),  and  g{z)  are  vector  valued  functions  and 


5(2)  =  fiz) 


(3.83) 


and  where  g{z)  satisfies  the  boundary  conditions  for  y(2).  Using  this  simplification  formula,  the 
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Green’s  function  for  the  fourth  moment  becomes 

G<4>  =  (^)"exp  -  B!)  •  (/  -  pO}  (3.84) 

/  /  D^(ri]D^(r2]  exp  <  jk  /  ri(T)  •  r2(T)<iT 

V  .'rj(0)=rj  L  Jo 

expj— f’[ri(T),r2(T),p(T)]dT| 

G<4>  is  seen  to  be  independent  of  the  value  of  p(t).  Because  of  this  independence,  the  final  value 
of  G<4>  will  not  change  irrespective  of  the  value  of  p(t).  Therefore  it  makes  sense  to  assign  it  a 
value  which  will  simplify  our  calculations.  For  reasons  which  will  soon  become  apparent,  we 
will  set  p(t)  equal  to  zero.  This  being  the  case,  it  follows  that  the  Green’s  function  becomes 

G<4>  =  (2^)^*?  •  (P"  -  /»')}  (3.85) 

f  ^  ’  D^[ri]D^[r2]exp  fjk  f  r,(T)  •  r2(T)dT 

J  .'rj(0)=rj  I  Jo 

expj— F(ri(T),r2(T)]dT| 

where  the  abbreviation  F[ri  (t),  r2(T)]  =  F’frj  (t),  r2(T), 0]  has  been  used.  At  this  point,  we 
recall  that  the  Green’s  function  was  composed  of  the  effects  from  4  independent  point  sources  or 
beams.  Now,  we  see  that  the  Green’s  function  has  been  manipulated  into  a  form  that  is  equivalent 
to  a  Green’s  function  with  two  beams.  Therefore,  the  Green’s  function  now  represents  two 
effective  point  sources  or  beams. 

Now  consider  the  decomposition  of  the  two  trajectories  rj  (t)  and  r2(T).  The 
decomposition  will  parallel  the  previous  decompositions  for  R{t)  and  p{t).  Specifically,  let 

»'i('r)  =  n('r)+  n  (t")  (3.86) 

r2ir)  =  rl{r)+r2{T) 


J  Jri(0)=r', 


r,(z)=i 


ir 

J  Jr\(0)=r'. 
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where 


r,(r)  =  (i^)r  +  ri 
»^(t)=  )-r  +  r2 


and  r'  =  »‘(t)|^_q,  t"  =  T'(T)|y_j,  and  so  fourth.  Additionally,  let 


n  (t)  =  £ 


^  ("r)  =  51  yn4>n{T) 


(3.87) 


(3.88) 


where  and  yn  are  the  ordered  pair  Fourier-sine  series  coefficients  and 


(t’’) 


(3.89) 


are  the  orthogonal  functions  on  [0,  z].  The  previously  mentioned  substitutions  for  ri  (t)  and  r2(T) 


will  result  in 


ri(T)  ■r2(T)  = 


(3.90) 


{(^)  +  (t")}  ■  {(^) + (tO) 

By  taking  the  dot  product  and  integrating  from  0  to  z  in  the  preceding  equation,  one  will  obtain  a 
simplification  of  Eq.  3.85  resulting  in 


(  ■Z  )  J  yr|(0)=r|  J  Jri(Q)=Tj 

{oo  Trik^  ^ 

J  51  •  yn  -  ^  F[tT(t)+  loo  (T),f2(T)+  iT),p{T)]dTj 


(3.91) 


where 


N 

(t)  =  Ylxn<i>niT)  (3.92) 

n=l 

N 

(t)  =  '^yn<f>niT) 

n=l 


and  where  i„,  i/„  6R^. 

Previously  we  noted  that  G<4>  is  independent  of  the  value  of  p(t),  and  we  choose  a 
convenient  value  for  it,  namely  zero.  We  now  wish  to  examine  the  consequence  of  such  a 
decision.  Recall 

+  (3.93) 

where 

p'  =  p\-f^2  +  Pi-PA  (3-94) 

J'  _  J' 

P  —  Pi  -  P2  P2  ~  Pa 

Let’s  consider  the  explicit  evaluation  of  the  the  preceding  two  equations  for  two  specific  values  of 
r.  When  r  =  0  we  have 


Pi  ~  P2  —  Pa  ~  Pi 


(3.95) 


and  when  t  =  zit  follows  that 


P'{  -  P2  =  PA  -  P'i  (3.96) 

Recall  from  elementary  vector  analysis,  two  vectors  (directed  line  segments)  are  equal  if  and  only 
if  they  are  parallel.  The  location  of  the  two  vectors  is  irrelevant.  Therefore,  the  choice  of 
P(t)  =  0  implies  that  the  vector  p\  -  P2  must  be  parallel  to  -  p'^  just  as  p"  —  pj  must  be  parallel 
to  /?4  -  p’{.  The  first  equality  states  that  a  quadrilateral  with  the  source  points  Pp  P3,  at 
each  of  its  vertices  is  in  fact  a  parallelogram.  The  same  holds  for  the  four  observation  points  ff{. 
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pl^,  f/^.  Fig.  3.1  illustrates  this  more  clearly.  This  geometry  of  source  and  observation  points 
will  not  create  any  problems  for  our  derivations  as  in  the  limit,  the  sides  of  the  parallelogram  will 
become  infinitesimally  small  and  the  four  poirv^;  will  be  merged  into  a  single  one  (24:56). 


It  would  appear  that  the  equation  for  the  fourth  statistical  moment  is  hopelessly 
complicated  with  the  infinite  number  of  integrations  contained  in  the  Feynman  integral. 

However,  as  in  the  case  of  the  single  path  integral,  we  will  approximate  the  infinite  number  of 
integrations  inherent  in  the  four-fold  two-dimensional  Feynman  integral  with  a  finite  number  (in 
reality  3)  of  ordinary  Riemann  integrals.  The  validity  of  this  approximation  may  be  verified 
either  by  mathematically  rigorous  means  (which  is  the  ideal  method)  or  by  comparing  the 
predictions  of  the  modified  equations  with  reality  (which  is  what  this  thesis  will  do).  As  our 
derivations  now  exist,  the  path  integral  //  D^[x]  is  the  integration  over  all  sinusoids  of  the  form 
given  by  Eq.  3.89.  The  approximation  of  the  path  integral  by  a  finite  number  of  functions  will 
result  in  a  finite  number  of  ordinary  Riemara;  integrations. 

As  the  interested  reader  will  find  in  Appendix  F,  it  is  possible  to  expand  the  functional 

/  F[r^{T)+x^{T)  +  r^{T),yoo{T)]dT  (3.97) 

Jo 

in  a  two-dimensional  variational  Taylor  series.  It  therefore  follows  from  Appendix  F  that  the 
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zeroth  variation  of  the  above  equation  for  N  terms  is  given  by 


f  F[riir)+Xff  (T),r2iT)+yN  {T)]dT  (3.98) 

Jo 

higher  order  variational  terms  will  be  defined  similar  to  the  one-dimensional  case.  As  already 
stated,  the  accuracy  of  our  approximation  will  be  validated  by  the  accuracy  with  which  the  results 
predict  re^ty.  Moreover,  the  task  of  actually  evaluating  the  resulting  integrations  must  be 
addressed.  As  it  was  stated  previously  and  as  it  will  be  repeated  here,  we  have  made  two 
approximations.  The  first  approximation  was  made  by  assuming  every  trajectory  of  the  Feyrunan 
integral  could  be  approximated  by  a  truncated  Fourier-sine  series.  The  second  approximation  was 
made  by  assuming  that  the  difference  between  the  7  functional  of  the  original  trajectory  and  the  7 
functional  of  the  approximate  trajectory  was  sufficiently  small  so  that  the  former  could  be 
approximated  by  the  later.  In  other  words,  the  approximation  by  the  first  term  in  the  variational 
Taylor  series  was  sufficiently  accurate.  The  case  of  JV  =  1  will  be  investigated  in  the  next  section. 

The  Fourth  Moment  for  a  Point  Source  In  this  section  we  are  now  ready  to  consider  the 
case  of  a  point  source  and  a  single  observation  point.  In  order  to  simplify  the  derivations,  we  will 
assume  that  both  the  source  point  and  observation  point  lie  along  the  z  axis  (tTC't)  =  r^ir)  =  0). 
Additionally,  we  will  consider  the  case  of  only  one  sinusoid  for  the  approximation  and  therefore 
the  symborwill  be  placed  over  a  variable  indicate  its  two-dimensional  nature.  Finally,  the 
turbulence  will  be  modeled  as  a  p-re  Kolmogorov  law.  Therefore,  the  fourth  statistical  moment 
apart  from  a  constant  becomes 

G<4>  =  K{z)JJ^^JJ^^Sxd^yexpi^jx  y-^J^  F[x<^(T),y<^(T)]dT|  (3.99) 

where  X  =  (x,,x2),y  =  {yuyi),<i>{T)  =  </>i(r)and 

F[x<f>{T),  y<f>{T)]  =  1.46C7f  {2\x\^/^  -I-  2(y|V3  _  ^(5/3  ^  ^  ^|5/3j  ^5/3^^)  (3  100) 
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and 


K{z)  =  ( «p  [(S"  -  Jf)  ■  (/>"  -  /)  +  M  -  r',)  •  (rj  -  r^)]} 


(3.101) 


Performing  the  integration  in  r  yields 


(3.102) 


G<4>  =  IC(z)  +X2y2)} 

(2(|i|=)>/‘  +  my'‘  -  (|i  -  6ff'‘  +  (\i  +  ift’'*)  } 


where  r(-)  is  the  Gamma  function.  The  previous  equation  was  written  in  a  form  to  facilitate  the 
substitution  of  polar  coordinates.  By  making  the  substitutions 


X  =  (xi,i2)  =  (riCos(di),ri  sin(fii))  (3.103) 

y  =  (yi,y2)  =  (r2COs(02),r2Sm(ff2)) 


Eq.  3.102  becomes 

K{z)G<4>=ff  dOi  f  f  de2  f  [  ridri  f  f  r2dr2exp{ir,r-2cos(0i  -  0iJi)lO4) 
J  Jo  J  Jo  J  Jo  J  Jo 

exp{C(2)  (2r;/^  +  2r5/3)} 

exp  |C(z)  ((rj  +  Tj  -  2rir2COs(^i  -  ^2)  -  +  2rir2COs(0i  -  ^2)))} 


where 


The  substitution  of 


C{z)  = 


-1.46v/2fc’/<»r2(4/3)  ^2.11/6 
r(8/3) 


<f>  =  $\  —  02 
</>'  =  02 


(3.105) 


(3.106) 
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will  permit  to  be  integrated  by  inspection,  resulting  in 


(jj.  \  2  r2ir  fOO  roo 

—j  2ir  I  d<p  I  ridrx  J  r2<ir2exp{jrir2COs(<^)}  (3.107) 

exp  {C(z)  (2rp  +  } 

exp  |C(2)  ^(rj  +  rl  -  2rir2  cos(^))*''®  —  (r?  +  rf  +  2rir2  cos(^))^/*)  } 


Further  simplifications  are  possible  by  considering  the  integration  of  <t>  on  [0, 7r/2]  instead  of 
[0, 27r]  and  recognizing  the  symmetry  of  G<4>  along  the  line  ri  =  r2.  These  substitutions  will 
now  be  performed  one  at  a  time.  First,  the  integration  range  will  be  divided  in  half.  Specifically, 
the  integrand  of  the  previous  triple  integral  may  be  considered  as  a  function  of  cos(<^)  only. 
Additionally,  the  previously  mentioned  integrand  is  an  even  function  in  cos(^).  Then  we  may 
write 


f  f{cos{(l>))d(l>  =  f  ficosi<l>))d<f>+  f  /(cos(^))# 

Jo  Jo  Jx 

=  /  f{cos{(l>))d(f>+  I  f(cosi(f>  +  ir))d<l> 

Jo  Jo 

=  /  ficosi<l>))d<f>+  f  f{-cos{<l>))d<l> 

Jo  Jo 

=  2  /  f{cos{<t>))d(f> 

Jo 


(3.108) 


A  similar  series  of  simplifications  may  be  made  which  will  result  in  the  integration  range  of  (f> 
being  reduced  by  half  again.  This  second  reduction  in  the  integration  region  will  cause  the  factor 
exp(jrir2  cos(^))  to  be  replaced  with  2cos(rir2  cos(^)).  Therefore,  it  follows  that 


fr/2  roo  rrl 

G<4>  =  87r/sr(2)  /  d<f>  I  r\dT\  I  r2<ir2COS  (rir2COs(^))  (3.109) 

Jo  Jo  Jo 

exp|C(2)  -  (rf  +  r-f  -  2rir2COs(</>))®/^  -  (r,  +  rf  +  2rir2COs(<^))^/®)} 


The  final  simplification  will  be  to  reduce  the  integration  region  0  <  ri  <  oo,  0  <  rj  <  oo  with 
0  <  r2  <  rj  <  oo.  Thus,  we  have  arrived  at  the  final  integral  form 
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rrll  foo  rrl 

G<4>  =  l6TrK{z)  I  d(f>  I  rjdri  /  r2rfr2COS  (rir2C0s(^))  (3.110) 

Jo  Jo  Jo 

exp  |C(z)  {2r\^^  +  —  {r\  +  rf  -  2rir2COs(^))^/®  -  {r\  +  rf  +  2rir2COs(^))^/*)} 

One  of  the  final  substitutions  which  may  be  made  concerns  the  function  C{z).  The  Rytov 
solution  for  the  relative  intensity  fluctuations  is  given  by  (25:103) 

Pi  =  Q.WIClk’l^z^^l^  (3.111) 

We  notice  that  the  units  of  C*  are  meters"*/^,  the  units  of  k  is  meters"*  and  z  is  simply  meters. 
Therefore,  if  we  make  the  substitution 

C{z)  =  -Q.6\90l  (3.112) 


the  analysis  of  Eq.  3.1 10  may  be  performed  under  dimensionless  conditions. 

Eq.  3.1 10  is  a  highly  oscillatory  integral  and  because  of  this  property,  special  techniques 
must  be  used  to  try  to  evaluate  it.  The  naive  method,  also  known  as  the  brute  force  method,  is  to 
simply  perform  the  triple  integral  indicated  above.  The  final  integral  over  r\  may  be  evaluated  by 
considering  a  sequence  of  partial  integrations  for  larger  and  larger  values  of  ri.  In  the  limit  as  ri 
tends  to  infinity,  it  is  hoped  that  the  partial  contributions  to  the  total  integral  tend  to  zero. 
However,  in  this  case,  hope  does  not  mirror  reality  and  for  a  finite  precision  machine  convergence 
is  nearly  impossible  to  attain.  In  order  to  accelerate  the  convergence  we  can  replace  the  integrand 
with  a  function  which  converges  to  zero  much  more  rapidly  (18:1).  Specifically,  we  can  subtract 
the  value  of  the  integrand  for  very  large  values  of  ri  and  then  add  it  back  again  later.  Wrth  an 
abuse  of  the  limr,-,oo  notation,  we  will  want  to  consider 

(l  —  lim  -I-  lim  \6'kK{z)  f  d<f>  (3.113) 

\  r,-»(»  rt^ooj  Jq 

roo  fT\ 

I  r\dr\  I  r2dr2COS  (rir2COs(<^)) 

Jo  Jo 

exp  |-0.619/3o  +  2r2'^^  -  (rj  +  »*2  ~  2rir2  cos(<^))*  -  (rf  +  rf  +  2rir2COs(^))«)  | 
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Therefore,  consider  the  case  where  limr,-n.oo>  This  equivalent  to  the  case  where  |x|  — »  oo  in 
Eq.  3.103.  Therefore,  as  ri  tends  to  infinity,  Eq.  3.1 10  is  fully  equivalent  to  Eq.  3.102  which  is 
shown  below  for  convenience. 


G<4>  =  K(z)  2  d^xd^yexp  {jixiyi  +  X2y2)]  (3.114) 

exp  {-0.619/?2  (2(|xp)V6  +  2(|yp)5/6  _  _  ^|2)5/6  +  (1^  ^  y|2)V6)} 

In  the  limit  as  x  tends  to  infinity  and  y  remains  unchanged,  the  integrand  of  Eq.  3.1 14 

exp  {-0.619y3^  (2(lxp)5/«  +  2(\y\‘^)^f^  -  (|x  -  +  (|x  +  (3.115) 


may  be  rewritten  as 

exp{j(xm+X2y2)-0.6l9f3l  (2(|xp)5/6  +  2(|y|2)V6  _  (|i|2)5/6  +  (|^|2)5/6)  J  (3  116) 

which  simplifies  to 

exp  {i(xiyi  +X2y2)  -0.619^^  (2|y|*^^)}  (3.117) 

Therefore,  we  will  be  interested  in  the  following  integration 

L{0o)  =  JJj^2  JJ^2  dxidyidx2dy2&^p  {i(xiyi  +  X2y2)  -  0.619/0^  (2|y|’/^)}  (3.118) 

Which  is  equivalent  to  the  integration 

L{PI)  =  dxidyidx2dy2Cxp  {j(xiyi  +  X2y2)  -  1.238/0^  (jy?  +  yfl^^®)}  (3.119) 

The  integrations  with  respect  to  X]  and  yi  may  be  performed  first  where 

f  exp(jiiyi)dxi  =  27r6(yi)  (3.120) 

J  —  OO 
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resulting  in 


L{01)  =  27r  dyi^Cl/i)  exp{ja:2y2}  (3.121) 

exp  I  - 1 .238/?^  (y?  +  yf )  |  = 

27r  J ^2  ‘ia:2dy2  exp{jX2y2}  exp  | - 1 .238j0o  (yf) | 

The  integrations  on  X2  and  y2  may  be  similarly  performed  yielding 

(3.122) 

Therefore,  it  follows  that  we  will  be  interested  in  the  evaluation  of 

rir/2  foo  fri 

\6TtK{z)  J  d<f)  J  r\dri  J  r2dr2  cos  (rir2 cos(^))  (3.123) 

[exp  |-0.619y3o  -  (r?  +  rf  -  2r\r2COs{<t>))^^^  -  (rf  +  rf  +  2rir2COs(^))*^®)}  - 

exp  {-  1.238/5o^r®/^}]  +  4jr^Ji'(z) 

The  preceding  integral  is  the  final  form  which  may  be  numerically  integrated,  although 
we  will  not  explicitly  do  so  as  we  could  not  obtain  a  convergent  solution  for  all  positive  values  of 
In  particular,  rapid  convergence  for  values  of  0^  >  10.0  was  obtained  using  Gauss-Kronrod 
integration.  The  corresponding  value  of  z  the  propagation  distance  is  given  by  simply  solving 

Q.ZOnClk’l^z^^l^  =  0^>  10.0  (3.124) 

which  results  in 

2  >  (32.57Cr2jb-V6)6/n  (3.125) 

Therefore,  given  a  particular  optical  wavelength,  and  a  value  of  Cc  =  C„  which  is  altitude 
dependent,  we  can  obtain  a  minimum  propagation  range  for  which  a  numerical  solution  to 
Eq.  3.1 10  will  rapidly  conveige. 

An  equation  for  the  scintillation  index  very  similar  to  Eq.  3.1 10  which  was  derived  under 
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a  different  set  of  assumptions  has  been  numerically  integrated  by  Beran  and  Whitman  (3:2140). 
Therefore,  the  result  of  performing  the  operations  of  Eq.  3.123  will  yield  a  function  of  Pq  which 
will  closely  match  the  scintillation  index  derived  by  Beran  and  Whitman.  Their  calculated 
scintillation  index  was  shown  to  be  accurate  to  the  theoretical  scintillation  index  within  ten 
percent.  Therefore,  we  expect  a  similar  magnitude  in  error  for  Eq.  3.123  in  its  predictions  of  the 
scintillation  index.  The  scintillation  index  from  Eq.  3.123  is  qualitatively  very  similar  to  the 
scintillation  index  for  strong  medium  fluctuations  which  appears  in  Fig.  2.4. 

This  chapter  has  examined  a  small  aspect  of  the  problem  of  wave  propagation  in  a 
randomly  inhomogeneous  medium.  The  next  chapter  will  examine  some  of  the  consequences  of 
this  examination. 


3-32 


IV.  Conclusions,  Contributions  and  Recommendations 


Research  Questions  Answered 

During  the  course  of  research  for  this  thesis,  it  was  desired  to  obtain  the  answers  to  some 
relevant  questions  which  are  repeated  from  Chapter  1  for  convenience. 

1.  Of  the  presently  known  subsets  which  have  any  practical  usefulness,  what  accuracy  do  they 
provide  when  used  to  help  approximate  functional  integration? 

2.  What  is  the  computational  complexity  of  these  subsets? 

3.  Do  the  characteristics  of  the  function  subsets  used  to  approximate  the  functional  integrand 
dictate  the  use  of  a  particular  numerical  integration  scheme? 

This  thesis  investigated  the  approximation  of  an  art>itrary  continuous  trajectory  with  a  truncated 
Fourier-sine  series  of  only  one  sinusoid;  therefore,  the  answers  to  the  above  three  questions  will 
have  to  be  accepted  with  such  an  understanding.  Specifically, 

1.  The  accuracy  of  a  single  sinusoid  is  estimated  to  be  on  the  order  of  ten  percent  error 
between  the  predicted  scintillation  index  and  the  theoretical  scintillation  index  for  a  point 
source  based  upon  the  results  of  Beran  and  Whitman  (3:2140).  This  estimation  of  the  error 
comes  from  the  comparison  of  the  three-fold  integral  obtained  from  the  Feynman  integral 
used  in  this  thesis  to  obtain  the  fourth  moment  of  the  field  and  the  similar  three-fold 
integral  examined  by  Beran  and  Whitman.  They  obtained  an  error  of  about  ten  percent  in 
their  research  and  its  reasonable  to  assume  the  error  in  this  thesis  is  on  the  same  order  of 
magnitude. 

2.  The  computational  complexity  of  this  particular  approximating  function  (a  single  sinusoid) 
was  extremely  high.  Initial  computational  runs  on  a  Cray  -YMP  using  Gauss-Kroiuod 
adaptive  integration  schemes  often  did  not  converge  in  under  one  hour.  Higher  values  of  /?q 
did  not  converge  at  all.  The  results  of  Beran  and  Whitman  were  obtained  by  heuristically 
choosing  particular  integration  regions  and  examining  the  limit  of  the  resulting  sequence  of 
partial  sums.  They  reported  acceptable  convergence  times  which  we  could  not  duplicate. 
However,  their  results  were  accepted  as  valid. 
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3.  Based  upon  the  experience  gained  in  writing  this  thesis  and  the  results  reported  by  Beran 
and  Whitman,  it  would  seem  that  the  analysis  of  the  integration  scheme  is  extremely 
important.  Apparently,  the  integration  region  must  be  subdivided  so  that  the  available 
integration  schemes,  for  example  Gauss-Kronrod,  will  achieve  acceptable  convergence  in  a 
reasonable  time.  The  exact  nature  of  the  subdivided  integration  regions  is  presently  very 
poorly  understood. 

Contributions 

Using  the  present  sinusoidal  approximation  to  the  path  integrand  and  the  presently 
implemented  numerical  integration  scheme  to  calculate  Eq.  3.1 10,  we  see  that  are  some 
Limitations  as  to  what  size  of  a  propagation  distance  will  yield  a  convergent  numerical  solution. 
Specifically,  as  of  the  time  of  the  writing  of  this  thesis,  the  limitation  is  given  by 

2  >  (4.1) 

meaning  that  values  of  z  larger  than  this  will  yield  a  convergent  solution.  As  an  example,  using  a 
measured  at  12497  meters  (above  sea  level)  of  5.0  10“*'*  (meter“^/^)  and  an  optical 
wavelength  of  7000  angstroms  (red  light),  we  find  a  convergent  solution  for  z  >  4550  meters 
(27:71). 

Recommendations 

In  order  to  obtain  a  better  understanding  of  the  method  of  path  integration  to  find  the 
scintillation  index,  it  is  recommended  that  the  nature  of  the  highly  oscillatory  integrand  used  in 
the  Eq.  3.123  be  better  understood.  Specifically,  the  necessary  integration  scheme  and  region 
over  which  such  a  scheme  is  applied  must  be  rigorously  understood.  Additionally,  a  single 
sinusoid  as  an  approximation  to  a  general  trajectory  is  clearly  not  a  good  approximation.  It  is 
possible  that  there  exist  other  orthonormal  function  sets  which  could  be  used  to  make  the 
necessary  approximations  (perhaps  wavelets).  Lastly,  though  this  woik  was  directed  at  gaining  a 
better  understanding  of  the  problem  from  the  path  integral  standpoint,  there  may  be  important 
future  applications  which  can  be  derived  from  the  work  contained  herein.  Therefore,  this 
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direction  of  research  should,  if  possible,  continue  to  receive  support 
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Appendix  A.  Parabolic  Wave  Equation 


As  a  convenience  to  the  reader,  this  appendix  has  been  provided.  This  appendix  will 
follow  the  method  used  in  (29:21 1-212)  for  the  development  of  the  parabolic  wave  equation. 
Starting  from  Maxwells’  equations  and  considering  a  non-conducting,  charge-free  medium  with  a 
unity  magnetic  permeability,  and  a  slowly  temporally  varying  dielectric  permeability 
e(r,  t)  =  c(p,  z).  Maxwell’s  equations  can  be  written 


V 


c  dt 


(A.1) 


V 


^  1  djeE) 
c  dt 


(A.2) 


V  •  =  0 


(A.3) 


where  E  is  the  electric  field  strength,  H  is  the  magnetic  field  strength,  r  =  (x,  y,  2)  is  a  position 
vector,  V  is  the  gradient  operator,  €  is  the  dielectric  permeability  of  the  medium,  and  c  is  the 
speed  of  light. 

In  the  atmosphere,  the  dielectric  permeability  is  very  slowly  changing,  therefore, 
e(r,  t)  ~  e(p,  2).  Thus,  letting  the  varying  portion  of  c  be  expressed  as. 


{^) 


(A.4) 


or  by  neglecting  the  explicit  (p,  2)  dependence 

7=lzJl> 

<«) 


(A.5) 


where  (•)  denotes  ensemble  averaging.  If  the  reasonable  requirement  that  the  atmosphere  exhibits 
non-relativistic  characteristics. 


{\v\) 


<  (I  «  l>, 


(A.6) 
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where  v  is  the  characteristic  velocity  of  atmospheric  turbulent  inhomogeneities,  and  c  is  the 
velocity  of  light  is  made  then  it  follows  that, 

(V^  +  €  E  -  V(-£^V  c)  (A.7) 

It  can  be  shown  that  the  last  term  in  the  above  equation  can  be  neglected  for  large-scale 
inhomogeneities:  Iq  This  being  the  case,  a  stochastic  wave  equation  which  holds  for  any 
component  of  the  electric  held  intensity  E  results. 


(V^  +  k^)E  =  -k^7  E 


(A.8) 


Thus, 


E  +  k\l+7)E  =  0 


(A.9) 


for  any  component  x,y,  or  z.  The  above  equation  is  known  as  the  stochastic  Helmholtz  equation. 
A  rigorous  derivation  of  the  parabolic  approximation  to  the  wave  equation  can  be  found  in 
(25:38-44).  However,  for  this  synopsis,  a  more  intuitive  and  simplistic  derivation  will  be  used. 
This  derivation  will  parallel  the  method  presented  in  (28:5.1 1).  An  optical  wave  propagating  in 
the  z  direction,  may  be  written  as. 


U  (i,  y,  z,  t)  =  u(a;,  y,  z)  cos(fcz  —  ujt) 


(A.  10) 


where  U (i,  y,  z,  t)  is  the  optical  field  at  the  point  (z,  y,  z)  and  the  time  t,  u{x,  y,  z)  is  the 
real-valued  amplitude  of  the  wave  at  the  point  (i,  y,  z),k  =  ^  and  A  is  the  carrier  wavelength, 
and  uj  is  the  angular  frequency  of  the  carrier.  By  considering  U  to  be  the  real  part  of  an  analytic 
signal,  Ua,  it  follows  that, 

I7(z,y,z)  =  ^{Ua{x,y,z)exp{jkz)exp{-jwt)}  (A.ll) 

where  S?(-)  =  \{i-)  +  (•)*).  and  (•)*  means  complex  conjugation.  By  considering  a 
monochromatic  optical  wave  and  the  variation  in  time  to  be  implicit,  it  is  possible  to  write  the  last 
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equation  as, 


(A.  12) 


U{x,y,z)  =  ^{Ua{x,y,z)exp(jk.z)}  =  ^{u{p,z)txp{jkz)} 

where  u{p,  z)  now  represents  the  complex  amplitude  of  the  optical  wave,  and  p  =  (a:,  y). 
Henceforth,  without  loss  of  generality,  this  thesis  will  consider  monochromatic  disturbances 
(waves)  which  propagate  in  approximately  the  z  direction.  Suppose  u(p,  z)  changes  little  on  the 
scale  of  A.  Then  it  follows  that  the  change  in  «(/>,  z),  which  is  represented  by  Su(p,  z),  will  be 
much  smaller  than  the  magnitude  of  u(p,  z).  In  other  words 

|<5u(/),z)|  <  !u(p,z)|  (A.13) 


Next,  consider  the  incremental  change  in  u(/>,  z)  along  an  incremental  change  in  z. 


oz 


(A.  14) 


Now,  because  it  has  been  assumed  that  u(p,  z)  varies  slowly  with  respect  to  the  carrier  frequency 
in  the  z  direction,  consider  £z  ss  A  in  which  case. 


du{p,z)  ^  ^  du(p,z) 

dz  dz 


(A.15) 


Thus,  it  follows  that. 


du{p,z)\ 
dz  I 


A  «:  \u{p,z)\ 


(A.  16) 


This  results  in. 


du{p,  z) 


dz 


y|«(p,z)|  <  k\u{p,z)\ 


Following  a  similar  line  of  reasoning,  a  similar  result  for  ^  follows. 


dz^ 


<  k 


dujp,  z) 
dz 


(A.17) 


(A.18) 


The  substitution  of  equations  A.  16,  A.17  and  A.18  into  equation  A.9  results  in  the  simplified 
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scalar-stochastic  parabolic  wave  equation, 


dz 


+  vi  +  « 


(p,  2)^  u(p,  2)  =  0. 


(A.  19) 


where  vl  =  ^  important  to  notice  that  equation  A.  19  is  a  first  order  partial 

differential  equation  in  z,  so  there  will  be  only  one  boundary  condition  in  z. 
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Appendix  B.  Method  of  Smooth  Perturbations 


Asymptotic  Expansion 

This  appendix  is  incuded  as  a  convenience  to  the  reader  who  may  not  be  fully  aware  of 
the  method  of  smooth  perturbations,  also  known  as  Rytov’s  method.  The  derivations  in  this 
appendix  will  closely  follow  those  found  in  (25:46-50). 

A  rigorous  derivation  of  the  parabolic  equation  would  exclude  the  effects  of 
backscattering  for  a  wave  propagating  in  an  inhomogeneous  medium.  Therefore,  Eq.  A.19  holds 
for  small  angles  about  the  mean  propagation  path.  The  problem  with  equation  A.19  is  that,  except 
in  the  most  trivial  of  cases,  it  is  not  solvable  in  closed  form.  Therefore,  the  method  of  smooth 
perturbations,  also  known  as  Rytov’s  method  is  applied.  This  method  will  approximate  the  actual 
solution  by  using  the  first  term  of  an  asymptotic  series  expansion.  Consider  an  equation  of  the 
form, 

/(x,O  =  0.  (B.l) 

where  ^  is  a  small-valued  parameter.  One  method  to  solve  the  above  equation  is  to  assume  a 
solution  of  the  form, 

fM)  =  'ZCMx)  (B.2) 

«=o 

where  {/i(a:)}“o  ^  sequence  of  (possibly  unbounded)  functions.  The  solution  to  B.2  is 

obtained  by  equating  similar  magnitudes  or  powers  of  For  example,  suppose  it  is  desired  to 
solve  the  equation 

exp(i^O  =  cos(i  4- 0  (B.3) 

for  small  values  of  By  inspection,  when  ^  is  zero  then  the  solution  is  simply,  1  =  0;  however, 
we  might  also  be  interested  in  solutions  in  x  for  small  non-zero  values  of  By  expanding  the 
exponential  and  sinusoidal  terms  into  their  respective  Maclurian  series  in  the  above  equation  we 
have 
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/  +  2xi 

V  ^ 


(B.4) 


21^  +  +  4x^^  + 


4! 


■) 


We  can  obtain  a  solution  by  equating  similar  powers  in  For  example,  the  first  few  simultaneous 
equations  will  be  obtained  by  equating  terms  in  and 


1  = 
x^  = 

2! 


x^  x'* 
'’2!  ■^4! 
2x  4x^ 

—  1  6x^ 


It  is  easy  to  see  that  the  first  of  the  preceding  three  equations  will  result  in  the  equality 


1  =  cos(x) 


(B.5) 


which  immediately  implies  x  =  2n7r  for  n  - - 2,  —  1, 0, 1, 2 •  •  •.  Each  equality  based  upon 

the  next  higher  power  in  ^  will  further  reduce  the  allowed  set  of  possible  solutions.  It  is  entirely 
possible  that  there  will  not  be  a  solution.  This  can  be  seen  to  occur  if  every  possible  solution  is 
eliminated  by  a  succeding  equation.  It  is  this  type  of  solution  technique  upon  which  the  Rytov 
method  is  based. 


The  Rytov  Method 

The  complex  amplitude,  u(p,  z)  of  a  propagating  wave  may  be  rewriten  in  polar  form. 


u{p,z)  =  >loexp 


,  Mp,z)  ,  c/  ^ 
In — - - hjS(p,z) 

^0 


=  Aq  exp(;f  +  jS)  =  Ao  exp[$(p,  z)]  (B.6) 


where  i4o  is  the  amplitude  of  an  initially  uniform  plane  wave  propagating  into  the  random 
inhomogeneous  medium  (25:48).  If  the  substitution  44ocxp[$(p,  z)]  =  u(p,  z)  is  made  into 
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equation  A.  19  one  will  obtain  the  equation; 


+  vi^  +  +  k^'^  {P, -e)  =  0  (B.7) 

This  equation  is  less  difficult  to  solve  than  equation  A.  19  as  the  fluctuating  term  e  as  an  additive 
term  and  not  a  multiplicitive  factor.  As  stated  in  (25:48),  equation  B.7  is  solved  by  assuming  an 
asymptotic  solution  of  the  form  $  =  $i  +  <&2  +  •  • ..  Additionally,  the  scale  size  of  is  assumed 
to  be  on  the  order  of  Thus,  |$i  |  >  |$2|  >  •  •  •  and  equation  B.7  can  be  recast  into  a 

collection  of  recursive,  though  linear  equations.  The  set  of  equations  was  generated  by  grouping 
terms  with  the  same  relative  magnitude.  Terms  which  are  on  the  size  of  a^,  and  so 

on,  are  collected  into  related  equations,  where  (Tt  is  the  standard  deviation  of  the  medium’s 
dielectric  permeability  e. 

+  €  (^,-z)  (B.8) 

+  Vi^2  *  -(Vx^i)^  (B.9) 

+  vi^3  =  -2  Vi  ^1  X  Vj.^2  (B.IO) 

Thus,  equation  B.8  is  of  size  a^,  equation  B.9  is  of  size  (cTe)^  and  equation  B.IO  is  of  size  (a^)^ 
(25:48).  What  usually  happens  in  the  method  of  smooth  perturbations  is  that  $  is  approximated 
by  only  the  first  term.  Thus,  $  ~  $i .  In  order  for  this  approximation  to  be  reasonable,  it  is 
necessary  for  I#2|  to  be  email.  This  will  happen  if  (vi(^i)^  is  small.  Equation  B.9  will  have  a 
solution  in  $2  of  the  form: 

fZ  fOO  fOO 

^2(p,z)=  /  /  (B.ll) 

Jo  J—00  J~oo 

where  the  Green’s  function  for  the  problem,  K  is  given  by 

If  ( Vi(^i is  small,  then  $2  in  equation  B.ll  will  be  small  and  the  approximation  $  ~  $1  will 


jk\p  -  p'\ 


/|2 


2iz-z') 


(B.12) 


dz 

d^2 

dz 


Ijk 

2jk 
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be  accurate.  The  magnitude  of  (Vj.(^i)^  necessary  for  the  above  approximation  to  hold  has 
been  found  to  be  (Vi(^i)^  <  If  iWs  holds  then  (Vx(^i)^  being  small  implies  transverse 
variations  in  $  must  be  relatively  small  and  therefore  $  must  be  relatively  smooth  in  the 
transverse  plane.  This  is  why  this  method  is  known  as  the  method  of  smooth  perturbations. 

If  Rytov's  method  is  applied  to  the  case  of  an  initially  unit  amplitude  plane  wave  striking 
a  phase  screen,  then  it  has  been  shown  that  the  phase  variance  of  the  resulting  wave  is  given  by 

I3l  =  0.3(nc^klz^i  (B.13) 

where  is  the  structure  constant  of  the  medium’s  dielectric  permeability,  fc  =  x  ^  I® 
carrier  wavelength,  z  is  the  propagation  distance  from  the  phase  screen  to  the  measurement  plane. 
One  should  immediately  notice  the  nearly  linear  dependence  in  k  and  nearly  quadratic 
dependence  in  z  for  the  parameter  I3q.  The  parameter  Pq  will  prove  to  be  very  useful  in  allowing 
the  path  integrals  to  be  expressed  in  a  dimensionless  form. 
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Appendix  C.  Method  of  Markov  Approximations 


This  appendix  has  been  provided  as  a  convenience  to  the  reader  and  draws  heavily  from 
(29:215-217).  It  will  provide  a  more  detailed  examination  of  the  Markov  approximation  method. 
Beginning  with  the  paraboLc  wave  equation  and  following  the  method  in  (25:85-86)  one  has 

(2jA:^  +  vi  +  «  {pyz))u{p,z)  =  0  (C.l) 


Now  we  must  recall  the  two  conditions  which  must  be  satisfied. 

First,  u(p,  z)  must  be  statistically  independent  of  the  effects  of  the  random 
inhomogenities  fo*"  any  location  z’  >  z  in  the  direction  of  propagation  (no  backscattering  effects). 

Second,  the  longitudional  correlation  radius  of  e  (p,  z))  must  be  less  than  the 
characteristic  time  of  changes  in  u(p,  z).  This  means  that  u{p,  z)  can’t  resolve  the  individual 
effects  of  the  fluctuations  in  c  {p,  z).  With  these  assumptions,  the  procedure  is  to  consider  two 
terms  of  equation  C.l  and  apply  the  identity: 


(^  +  /(p,2))v(p  «)  =  exp(-J^  (C.2) 


Therefore, 


(27*^  +  k^'e  (p,z)^  v(p,z)  = 


2jk 


{mf '  I  [“f  [wk  r ' 


(C.3) 


Thus,  we  have: 


2jk  |exp  c  (/),  ^  [exp  e  (/?,  v{p,  z)  =  -  Vx  ‘«ip,  z) 


(C.4) 


Rearranging  the  factors  yields 


2jk-^  exp(jf  (p,^)d^v{p,z)  =  exp^yjT  e  vl v(p,z)  (C.5) 


C-1 


Next,  integrating  from  0  to  z  yields 


Continuing, 


-2jk  |exp  e  v(p,z)  -  2jkv(p,0)  = 

exp  c  (P,0<^)  «(P,^)  = 


(C.6) 


(C.7) 


Solving  for  v{p,  z)  yields. 


v{p,z 

'  =  [ 

-  exp  y  «  v(p,0)  + 

.1“'’ 

(-y  vi  v{p,z')dz' 

(C.8) 


Now,  the  last  two  integrals  in  the  above  equation  may  be  combined  to  yield, 
vip,z)  =  |-exp  ^  v(/7,0)+ 


(C.9) 


Now,  the  expectation  of  both  sides  is  taken. 


{vip,z))=  (C.IO) 

iwk  r®”  (/' 


Consider  the  second  term  in  the  last  equation.  Its  expectation  may  be  simplified  by  noting  that  the 
factor  //,  e  {p,Od^  is  dependant  only  on  values  of  z  for  z  >  z’,  while  the  factor  viw(/J,  z')  is 
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dependent  only  on  values  of  z  for  z  <  z’.  Therefore,  the  two  factors  are  statistically  independent 
except  on  a  set  of  measure  zero.  So,  the  expectation  operation  may  be  simplified. 


<«(/>,  «)> 

-1 

23k 


r  (r  f  ‘  <'>■«'«))  (vi-ta 


(C.11) 


If  the  statistics  of  e  (p,  0  ^  known  then,  equation  C.l  1  is  a  closed  equation  for  the  mean  of  the 
field  amplitude.  For  example,  the  statistics  of  ?  (p,  0  might  be  gaussian.  If  this  is  the  case,  a 
closed  partial  differential  equation  for  the  mean  of  v{p,  z)  can  be  obtained  (25:87).  Higher  order 
moments  of  v(p,  z)  may  be  obtained  by  generalizating  the  above  method  (25:88)  and  this  will  be 
done  shortly. 

By  performing  steps  similar  to  those  already  demonstrated  it  is  possible  to  derive  a  closed 
equation  for  an  arbitrary  statistical  moment.  In  particular  if 

rn,m('2)  —  rn,m(Pl  ifKli'"  i  Pni  P\i  PTi't ' ' '  >  Pm^^)  ~  (C.12) 

(u(pi,  2)u(P2,  2)  •  •  •  «(p„,  2)u*(pi,  2)«*(p2,  2)  •  •  •  Z)) 


and 


—  Qn,m (Pl  ?  P2i  '  '  '  i  Pni  PliP2i '  '  '  i  Pmi  — 
e  e  (P2,2)  +  ---+  €  iPn,z)-T  {p\,z)-T  (p^,2) 


(C.13) 


then  as  shown  in  (29:216),  it  follows  that 

2jkTn,m{z)  -  2jkm  ^exp  Q(x)di|^  +  (C.14) 

Jo  Jo  (vi  +  •  •  ■  +  Vn  -  V? - Vm)  r„,,„(l)dx  =  0 

where  V?  =  ^-  This  is  a  closed  partial  differential  equation  for  an  arbitrary  statistical  moment 
of  the  complex  amplitude.  If  a  further  assumption  is  made  about  e ,  namely  that  e  has  a  Gaussian 
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distribution  and  is  delta  correlated,  specifically 


G  ipi .  )  «  (p2, 22))  =  A(p\  -  Pi)6G\  -  ^2) 


(C.15) 


then  as  shown  in  (29:217)  it  follows  that 


2jk^r„,m(z)  +  (vi  +  Vn-V? - Vm  •  •  •)  Tn.mC^)  +  ^i=’n,mr„,m(z)  =  0 

(C.16) 

where 


7^n,rn  —  n,m(Pli  P2f  '  '  '  i  Pni  P\i  P2i  “  '  t  Pm)  — 

n  n  n  m  mm 

51 51  -  p6)  -  2  5Z  5Z  -  pfc)  +  51 -  p'h) 

a=:lb=l  a=l  b=l  a=16=l 


(C.17) 


The  assumption  that  c  is  Gaussian  distributed  follows  from  the  law  of  large  numbers  which  states 
that  the  sum  of  independent  and  identically  distributed  random  variables  tends  to  a  Gaussian 
distribution.  The  assumption  of  delta  correlation  is  a  good  approximation  if  the  correlation  radius 
of  "t  is  small.  The  solution  for  r4(z,  pi,P2,p3,P4)  can  be  obtained  for  a  phase  screen  of  thickness 
(^z)  and  is  found  to  be 


(C.18) 


exp  {^l(Pl  -  p\  P  -  |p2  -  Plf  +  l/>3  -  P3p  -  \p4  -  Pa?^ 

exp  I  ^  (2D(p',  -  p'2)  -  lD{p\  -  p\)  +  lD{p\  -  p',)  +  2D(p^  -  p^))  | 

exp  {2D{p'j  -  pi)  +  2£)(p'  -  p^))  |  d^p'.d^p'^d^p'^d^p'. 


where  D(x)  =  >1(0)  -  >l(x).  For  an  initially  totally  coherent  plane  wave,  r4(pj ,  p2,  P3,  P4, 0)  = 
1.  A  final  simplification  may  be  perfomed  by  using  the  substitutions 


p'l  -P2  =  n 

P\  -  Pa -'^2 


(C.19) 
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Pi -(^  =  *•3 
Pi +P^  =  2r 


Pi-  P2=  Pi 
Pi  -P4=  P2 


and  integrating  over  r  and  r3  results  in 

r4(«,Pi,P2)  =  (2^)  j  j &xp^Y^ri-pi)  {r2-p2)-^^F{ri,r2)^<fri<fr2 

(C.20) 

where  /’(ri,r2)  =  2D(ri)  +  2D{r2)  -  2D{t\  -  T2)  -  2D{r\  +  r2).  The  preceding  result  gives 
the  fourth  statistical  moment  at  a  distance  z  behind  a  phase  screen  of  thickness  6z  for  an  initially 
totally  coherent  plane  wave.  It  is  possible  to  consider  Eq.  C.20  another  way.  Examination  of 
Eq.  C.20  reveals  the  integration  of  a  function  over  the  entire  plane.  The  integrand  of  Eq.  C.20 
may  therefore  be  considered  to  be,  apart  from  a  constant,  a  (ireen’s  function  which  describes  the 
fourth  statistical  moment  of  a  point  source.  This  Green’s  function  has  the  form 

G'<4>(Pl,P2.P3t/>4)  =  (C.21) 

(7  As  k'^Sz 

^—{Pi  -  P%)-{p2-  Pa) - ^  [2D(pi)  +  2D{p2)  -  D{pi  -  P2)  -  D{pi  +  P2)] 

It  will  be  found  to  be  useful  to  compare  the  preceding  Green’s  function  with  a  functionally  similar 
path  integral.  This  is  done  in  chapter  3. 
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Appendix  D.  The  Feynman  Integral  Schrodinger  wave  Equation  Relationship 


Introduction 

This  appendix,  which  draws  directly  from  (11 :26-78)  is  included  to  help  the  reader 
understand  the  origin  of  the  Feynman  integral.  This  appendix  will  closely  follow  the  original 
method  of  deriving  the  Schrodinger  wave  equation  from  its  Feynman  integral  solution.  There  are 
several  reasons  for  including  this  appendix,  namely: 

1.  To  provide  the  reader  of  this  thesis  with  the  methodology  used  in  some  of  the  quantum 
mechanics  related  literature.  This  is  justifiable  because  as  previously  noted,  many  of  the 
methods  and  solutions  of  quantum  mechanics  have  been  used  to  solve  wave  propagation 
problems. 

2.  To  provide  an  alternative  example  of  the  relationship  between  a  parabolic  equation,  the 
Schrbdinger  wave  equation,  (a  first  order  in  time,  second  order  in  space  partial  differential 
equation)  and  a  Feynman  integral. 

3.  To  provide  some  additional  justification  to  use  the  functional  integral  approach  to  solve 
related  problems. 

The  Classical  Trajectory 

In  the  everyday  experiences  of  people,  the  effects  of  the  Hcisenbeig  uncertainty  principle 
can  largely  be  ignored.  That  is,  the  kinematics  of  common  objects,  cars,  horses,  books,  all  obey 
Newton’s  laws  of  motion.  The  expected  path  or  trajectory  that  a  moving  object  will  follow  is 
called  the  classical  path  (11 :26).  In  the  quantum  mechanical  world,  the  Heisenberg  uncertainty 
principle  states  that  possibilities  other  than  the  classical  or  expected  trajectory  are  possible.  Fig. 
shows  some  of  these  paths.  So.  according  to  the  Heisenberg  uncertainty  principle,  it  is  possible 
for  a  car  spontaneously  decay  or  for  a  horse  to  fly  to  the  moon.  Of  course  this  is  obviously  not 
how  we  have  come  to  expect  cars  and  horses  to  behave  and  therefore  the  probability  of  such 
occurances  is  vanishingly  small.  Thus,  the  classical  trajectory  is  in  the  quantum  mechanical 
world,  simply  the  trajectory  with  the  greatest  probability  of  occuring.  It  turns  out  that  there  is  a 
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functional  defined  on  all  the  possible  trajectories  which  exhibits  an  extremum  for  the  classical 
trajectory.  This  is  the  Action  functional. 


-  Possible  trajectory 

Classical  trajectory 

Figure  D.l.  The  subdivision  of  a  trajectory 


The  Principle  of  Least  Action 

Consider  a  particle  initially  at  (xa,  j/a>  ^a)  at  the  time  ta  which  travels,  in  a  continuous 
fashion,  to  (xb,  Vb,  ^b)  at  the  time  4  and  obeys  Newton’s  laws  of  motion.  The  continuous  path  is 
therefore  be  described  as  the  continuous  parametric  function  (x(t),  y{t),  z{t)).  The  action 
functional  of  the  particle  along  the  path  (x(^),  y(t),  z{t))  is  defined  as. 


5[x,  y,  2]  =  [Ek  -  Ep)  dt  =  I'"  L  x(t),  y(t),  2(t),  t)  dt  (D.  1) 


where  5[x,  y,  z]  is  the  action  fuiKtional  along  the  path  (x(t),  y(t),  z(t)). 


Ek  =  f 


(my 


+ 


the  kinetic  energy,  Ep  =  V (x(t),  y(t),  z{i),  t)  is  the 
potential  energy,  L(^^,  ^^,x(t),y(t),  z(t),  t)  is  the  difference  between  the  Newtonian 

kinetic  and  potential  energies.  The  action  functional  for  different  trajectories  can  either  be  similar 
or  different.  In  the  case  of  the  classical  trajectory,  the  action  functional  will  be  an  extremum.  This 
means  that  for  a  slight  perturbation,  6  of  the  classical  path,  there  should  be  no  difference  in  the 
action  to  the  first  order  in  S.  Thus,  if  (^x.  Sy,  Sz)  is  a  slight  perturbation  of  the  classical  path 
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which  matches  the  classical  path  at  the  end  points,  as  shown  in  Figthen  this  implies. 


5[x(t)  +  Sxit),y{t)  +  Sy{t),  z(t)  +  6z{t)\  -  5[x(t),  y(t),  z{t)]  =  0  (D.2) 


Classical  path 
Path  perturbation 


Figure  D.2.  A  path  with  and  without  a  slight  perturbation 


It  is  possible  to  expand  5[x(t)  +  6x{t),  y{t)  +  Sy{t),  z{t)  +  Sz{t)]  in  a  Taylor  series  about 
(x(t),  y(t),  z(t))  which  means  that  L  can  be  expanded  to  first  order  as. 


dx{t)  ,  ^dxjt)  dyjt)  .  ^dyjt)  dz{t)  .  ^dzjt) 
dt  dt  ’  dt  dt  '  dt  dt  ' 
x{t)  +  Sx(t),y(t)  +  Sy(t),  z(t)  +  Sz(t),  t) 

=  L(x  -hSx,y  +  Sy,  z  +  Sz,x  +  6x,y  +  6y,  z  +  dz,  t) 


ci:  L{x,y,z,x,y,z)  + 

dL^..dL^.  dL^..dL  dL 
Sx—  +  Sx—  +  6y—  +  6y—  +  6z—  +  6z— 
dx  dx  dy  ay  dz  dz 


(D.3) 


where  the  over-dot  represents  partial  differentiation  with  respect  to  t, second  and  higher  order 
terms  in  and  have  been  dropped.  Now,  the  action  functional  of  the  perturbed  path  is 
written. 


5[x  -I-  5x,  1/  -f  6y,  z  +  8z\  = 


'Sx^  +  6x^) 
<  ax  dx / 

f^.dL  ,  dL\ 


dt  -f- 
dt 


(D.4) 
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it  then  follows  that. 


SS  =  5[i  +  6x,y  +  Sy,  z-\-6z]-  5[a:,  j/,  z]  =  (D.5) 

Additionally,  it  is  necessary  that  8x  =  8y  =  8z  =  Oat  the  trajectory  endpoints  <„  and  tb.  Thus,  in 
order  that  £5  =  0  to  first  order  in  8x,  8y,  8z,  it  is  sufficient  that, 

/d  /dL\  dL\  fd  fdL\  dL\  (d  {dL\  dL\  ^ 

VdtVSiJ  dx  J  \dt\dyj  dy  J  \.dtv5i/  dz  J 


which  is  seen  to  be  the  classical  lagrangian  of  the  Newtonian  particle.  The  notion  of  a  trajectory 
or  path  can  be  more  generalized  as  it  will  be  shown  in  a  later  section. 


Wave  Functions 

Nothing  in  the  quantum  mechanical  world,  as  already  noted,  is  certain.  Therefore  highly 
probable,  probable,  possible,  improbable  and  highly  improbable  events  all  have  associated 
probabilities.  A  wave  function  is  a  function  which  describes  these  probabilities.  Written  as 
V'(x,  t),  its  magnitude  squared  gives  the  probability  of  a  certain  event  at  the  location  x  at  the  time 
t.  That  is,  if  ip{x,  t)  is  the  wave  function  associated  with  an  event  A,  then  the  probability  of  event 
A  is 

(D.7) 

here  x  is  the  one-dimensional  space  variable.  It  is  easy  to  extend  x  into  three  dimensions 
(x,  y,  z).  Given  an  event,  there  are  several  possible  paths  or  trajectories.  Each  possible  path  or 
trajectory  has  an  associated  probability  of  occurance  which  in  turn  is  described  by  a  wave 
function.  If  there  are  two  sequential  events,  event  A  follower  oy  event  B,  then  there  is  composite 
wave  function  to  describe  the  composite  event  AB.  Furthermore,  there  is  an  important  relation 
between  the  individual  wave  functions  for  the  individual  events  A  and  B  and  given  by, 

roc 

1p(Xb,ti]  =  K{XbJb^Xa,ta)lp{Xa,ta)dXa  (D.8) 

7  —  00 
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where  K{xi„  tb,  Xa,  to)  is  another  wave  function  which  will  be  described  in  the  next  section. 


The  Wave  Function  K 

The  wave  function  K{x2,t2,xi,ti),  which  appeared  in  the  previous  section  is  an  integral 
over  paths.  It  is  given  by, 


I'’-  “I = !“„  i  /  /  ■  ■  ■  /  “"’  d®’''’  “0  ^  X  ^ ' 


dx„- 


n— 1 


(D.9) 


where  A  =  and  e  and  Xi  are  shown  in  the  following  figure,  and  the  dependence  of  K  on 

ta  and  4  is  understood. 


Figure  D.3.  The  subdivision  of  a  trajectory 


Often  one  finds  in  the  current  literature  the  approximation  of  continuous  trajectories  or  paths  by 
piecewise  linear  approximations  as  shown  in  Fi^29:226).  These  approximating  functions,  which 
are  also  called  polygonal  functions,  are  used  to  simplify  the  limiting  process  in  Eq.  D.9.  It  is 
important  to  note  that  given  any  continuous  curve  in  R^,  it  is  possible  to  approximate  it  as  closely 
as  desired  with  a  polygonal  function.  That  is,  given  a  continuous  curve  in  R^,  there  is  a  sequence 
of  polygonal  functions  which  converges  to  the  given  continuous  curve  and  infact,  the 
convergence  is  uniform.  We  now  have  all  the  necessary  tools  to  derive  the  Schrddinger  wave 
equation.  This  will  be  done  ir  the  following  section. 
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Figure  D.4.  The  approximation  of  a  trajectory  by  a  polygonal  function 


Derivation  of  the  Schrodinger  Wave  Equation 

The  one  dimensional  version  of  the  Schrodinger  wave  equation  will  now  be  derived.  As 
previously  discussed,  the  derivation  of  this  equation  will  help  reinforce  the  readers  understanding 
of  the  relationship  between  a  partial  differential  equation  and  its  path  integral  solution.  Beginning 
with  Eq.  D.8,  for  a  Newtonian  particle  at  position  Xa  at  time  ta  and  position  xi,  at  time  ta,  where 
Xa  and  Xb  are  points  on  the  continuous  trajectory  i.  Let  Itt  -ta\  =  €'  which  is  very  small  then, 


fOO 

J  —  00 
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For  values  of  e'  small  enough,  that  is,  smaller  than  e,  the  path  integral  in  Eq.  D.9  will 
reduce  to  a  regular  integral, 

K{xb,Xa)  ~  ^y^exp  (^e'^S[xb,Xa]^  dxa  (D  ll) 


For  small  t'  the  time  derivitive  of  the  continuous  trajectory  x  with  respect  to  t  may  be  written. 


X 


lim 

('—0 


fXb-Xa\  (Xb-Xa\  f  Xb  -  Xa 


) 


(D.I2) 


Additionally,  for  small  values  of  e',  Xb  is  close  to  Xo,  and  therefore,  with  error  on  the  order  of  e' 
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that  is  ©(O,  xb  c-  ~  ^^,t).  Therefore, 


or  1  ^  r  f  ■♦N 

—  i  I  ,  2 


(D.13) 


Therefore,  the  substitution  of  Eq.  D.13,  Eq.  D.l  1  and  lx/,  —  Xal=  t  into  Eq.  D.IO,  will  result  in 
Eq.  D.14, 

V'(x6,t  +  e')  =  y  ^  exp  e'^L  j  t/)(xa  +  t)dxa  (D.14) 


The  Lagrangian  for  a  Newtonian  particle  traveling  on  a  one  dimensional  trajectory  is 
given  by  Eq.  D.l 5 

.  (D.15) 


L(x,  X,  t)  =  -j-  +  V  (x,  t) 


Thus,  making  this  substitution  into  Eq.  D.14  results  in  Eq.  D.16, 


/oo  I 


exp 


,1  I  mr^  .  T  , 


1p(Xa,t)dXa  (D.16) 


This  may  be  simplified  to. 


V’(l6,f + 


1p{Xa,t)dXa 


(D.17) 


The  first  exponential  term  in  the  integrand  in  Eq.  D.17  is  highly  oscillatory.  Moreover, 
the  contrabutions  will  be  vanishingly  small  for  large  values  of  t.  Consider  values  of  t  such  that 
varies  by  at  most  one  radian  then, 


jmr'^ 

'me' 


<  I 


(D.18) 


Therefore,  values  of  t  such  that  Eq.  D.18  holds  will  give  a  non-zero  contribution  to  Eq.  D.17. 
A  function  f{x,  t)  which  is  continuous  in  all  differential  orders  may  be  expanded  in  a 
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Taylor  series  about  the  temporal  point  t'  as. 


^  [(a?) 


(D.19) 


By  considering  the  Taylor  expansion  of  the  wave  function  ipixb,  t  +  e')  about  the  point  t,  it  may 
be  approximated  to  first  order  in  e'  as. 


i/)(xb,  t)  ~  tpixf,,  t)  +  (D.20) 

Similarly,  exp  {xa  +  f ,  <))  may  be  expanded  in  a  Taylor  series  to  first  order  in  c'  as, 

exp  d  +  ^ViXa  +  ^,t)  (D.21) 

Lastly,  ip{xb,  t)  may  be  expanded  in  a  Taylor  series  about  Xa  to  first  order  in  d  which  by  Eq.  D.18 
means  to  second  order  in  t  as. 


ll>{Xb,t)  V’(a!a,t)  +  T— 


(D.22) 


X\,-Xa 


The  previous  approximations  may  now  be  used  in  Eq.  D.17  to  yield. 


04) 

/. 


■4>{xb,t)  = 


(D.23) 


1  /  irriT^ 


l-^-iv{Xa,t) 


d  Sd 

1  +r- - h 


dxb  2  dxl 


6J 


xl}{Xb,t)dT 


This  equation  can  be  simplified  by  noting  that. 


and  that. 


(D.24) 


(D.25) 
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By  equating  terms  in  Eq.  D.23  of  order  e'®  Eq.  D.26  results. 


This  will  hold  if  ^  By  using  the  simplifications  affored  by  Eq.  D.24  and  Eq.  D.25 

with  the  value  for  A,  the  differential  equation  of  Eq.  D.27  results. 


(i  + 1)  =  (i  -  (i  + 1^^)  (d.27) 


By  using  some  simple  algebra  and  remembering  that  e'r^  is  of  order  c'^  and  can  be  dropped, 
Schrodingers  wave  equation  results, 


dt 


(D.28) 


where  Xb  has  been  replaced  with  x.  Examination  of  the  above  partial  differential  equation  reveals 
that  it  is  a  first  order  in  time,  second  order  in  space  equation.  Although  it  is  not  elementary,  it  is 
possible  to  follow  the  above  steps  by  considering  the  space  variable  to  be  two  dimensional.  In 
that  case,  it  will  be  possible  to  demonstrate  a  relationship  between  the  two  dimensional 
Schrodinger  wave  equation  and  its  two  dimensional  Feynman  integral  solution.  Lastly,  by 
allowing  the  temporal  variable  to  be  replaced  with  a  one  dimensional  space  variable  which  is 
transverse  to  the  other  two  space  variables,  it  is  possible  to  convert  the  two  dimensional 
Schrodinger  wave  equation  into  an  equivalent  parabolic  wave  equation  which  will  have  an 
associated  Feynman  integral  solution. 
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Appendix  E.  Solution  of  the  Helmholtz  Equation 

This  appendix  is  included  as  an  aide  to  the  reader  who  may  wish  to  briefly  review  the  Helmholtz 
equation  and  some  of  its  characteristics. 

The  Solution  to  the  Deterministic  Helmholtz  Equation 

This  section  will  give  a  solution  of  the  deterministic  Helmholtz  equation  under  the  condition  of  a 
slowly  varying  complex  amplitude  in  the  z  direction.  That  is,  the  assumptions  for  the  parabolic 
approximation  hold.  In  this  case,  the  solution  to  the  Helmholtz  equation  Eq.  E.1  with  a  unity 
boundary  condition, 

{^^  +  k^)u{p,z)  =  0  (E.1) 

where  u{p,  z)  is  the  complex  amplitude  of  the  optical  wave  after  a  propagation  distance  of  z,  k  is 
^  and  A  is  the  carrier  wavelength,  p  =  (i,  y),  namely  p  is  in  the  plane  transverse  to  the  mean 
direction  of  propagation,  z.  A  solution  is  given  by  a  superposition  integral  of  the  form, 

uip,z)  G(p,/,z)u(p',0)dV  (E.2) 

where  G{p,  p',  z)  is  a  Greens  fiuicron  which  satisfies  the  equation 

{V^-^-k^)G{p,p',z)  =  6{p-p')  (E.3) 


where  6{p  —  p')  is  the  Dirac  delta  distribution  is  given  by  (4:378),  (25:56), 


d_  / exp(j7:^|p-/p-f  zA 


'|p-  p'P  +  z^ 


where  u(p,  z)  is  the  complex  amplitude  of  the  optical  wave  after  a  propagation  distance  of  z,  k  is 
^  and  A  is  the  carrier  wavelength,  p  =  (x,  y),  namely  p  is  in  the  measurement  or  observation 
plane  transverse  to  the  mean  direction  of  propagation,  z,  p'  =  (x',  y'),  namely  p  is  in  the  object 
plane,  z  =  0,  is  the  Greens  function,  E  is  the  region  where  u(p',  0)  is 

non-vanishing. 
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If  Eq.  E.4  is  considered  for  values  of  z  such  that  z  »  A  and  z»\p-  then  the  solution  to  the 
Helmholtz  equation  Eq.  E.l  may  be  approximated  by, 

(B.) 

Eq.  E.5  is  known  as  the  Huygens-Fresnel  integral. 

Interestingly  enough,  Eq.  E.5  is  the  solution  to  the  parabolic  approximation  to  the  deterministic 
Helmholtz  equation  Eq.  E.6, 

{2jk^^+v\)u{p,z)  =  Q  (E.6) 

where  u(p,  2)  is  the  complex  amplitude  of  the  optical  wave  after  a  propagation  distance  of  z, 
vi  =  ^  i  Thus,  the  solution  to  the  parabolic  equation  requires  the 

consideration  of  small  angles  about  the  mean  direction  of  propagation  which  is  the  same  as 
requiring  that  u{p,  z)  vary  slowly  in  the  mean  propagation  direction,  z.  But,  this  simply  means 
that  u(p,  2)  has  a  narrow  spatial  spectrum  for  the  z  spatial  frequency.  The  effects  of  the  medium 
have  been  implicitly  ignored.  However,  the  inclusion  of  the  random  process  e(p,  z)  will  not 
appreciably  complicate  the  solution  to  the  Helmholtz  equation  previously  derive. 

The  Solution  to  the  Stochastic  Helmholtz  Equation 

The  stochastic  Helmholtz  equation  has  already  been  derived  in  Appendix  A  from  Maxwell’s 
electromagnetic  field  equations.  It  was  found  to  be, 


(v'  +  fc"(l+ e  (^,z)))u(p,z)  =  0  (E.7) 

If  the  solution  to  Eq.  E.7  is  considered  for  the  region  such  that  the  Fresnel  approximation  holds 
then,  the  solution  will  take  the  form, 

«(/>,  z)  =  J  j u(p,  z)  exp  exp  e  {p' ,  z)^  d^p'  (E.8) 
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where  the  approximation, 


(E.9) 


has  been  made.  Additionally,  in  this  approximation  is  the  implicit  assumption  that  e  {p,  z)  is 
relatively  constant  over  a  small  increment  in  the  z  direction.  Thus,  the  above  approximation  will 
become  more  accurate  as  the  incremental  change  in  z  is  small. 
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Appendix  F.  Variational  methods 


Functions  and  Functionals 

This  appendix  which  closely  follows  (24)  is  included  as  an  aide  to  the  reader  who  may 
not  be  fully  conversant  in  functional  or  variational  differentiation.  The  terms  rule,  function, 
mapping  are  all  synonymous  in  this  appendix,  as  are  the  terms  variational  differentiation  and 
functional  differentiation.  In  addition,  a  great  amount  of  mathematical  rigor  has  been  sacrificed 
for  both  increased  readability  and  decreased  length.  If  this  loss  of  rigor  causes  excessive  concern, 
then  you  are  advised  to  refer  to  a  good  book  on  functional  and  variational  principles. 

Recall  that  the  simplistic  idea  of  a  function  is  a  mapping  or  rule  from  one  set  to  another. 
For  example,  the  function  given  by  y  =  2a:  looks  iimocent  enough.  However,  if  we  assume  x  is 
only  defined  on  the  set  of  integers  {1,2, 33, 44},  then  y  may  only  take  on  the  values  {2, 4, 66, 88}. 
Therefore  the  domain  of  the  function  (the  allowed  values  of  i)  and  the  range  of  the  function  (the 
allowed  values  of  y)  are  also  extremely  important  in  the  total  idea  of  the  function.  Taking  the 
simplification  of  the  idea  of  a  function  still  further,  we  could  consider  a  function  to  be  a  machine. 
We  send  numbers  into  the  machine  and  other  numbers  come  out.  The  types  of  numbers  we  can 
send  into  the  machine  is  constrained  by  the  function’s  domain.  The  numbers  the  machine  will 
return  are  constrained  by  both  the  domain  and  the  range.  Similar  to  the  idea  a  function  is  the  idea 
of  a  functional.  A  functional  is  a  rule  or  mapping  that  is  defined  on  some  set  of  functions  (the 
domain)  and  which  returns  real  numbers  (the  range).  A  functional  might  be  defined  by  the  rule 
which  returns  1  for  all  continuous  functions  and  0  for  all  discontinuous  functions.  Another 
example  of  a  functional  is  the  area  under  a  continuous  curve.  This  similarity  between  functions 
and  functionals  has  some  interesting  implications.  Therefore,  a  functional  is  also  a  machine; 
however,  the  functional  machine  has  functions  as  its  input  and  real  numbers  as  its  output. 


Function  Differentiation 


The  derivative  of  a  function  /  at  a  point  x  is  given  by 

dfjx)  _  fix  +  Ax)-f{x) 
dx  Ax— >0  Ax 


(F.l) 


whenever  the  limit  exists.  For  example  the  derivative  of  at  the  point  1  is  given  by 

Um  (F.2) 

dx  Ax-*o  Ax 

X=1 

Expanding,  cancelling,  and  evaluating  the  resulting  terms  results  in  the  expected  result  3. 
An  example  of  a  function  without  a  derivative  is  the  mapping  is  the  general  polygonal  function 
which  is  used  to  define  a  trajectory  in  its  limit  (see  appendix  D).  On  the  trajectory  where  two 
straight  lines  join,  the  function  describing  the  trajectory  will  not  be  differentiable.  The  single 
exception  is  the  case  when  both  lines  have  the  same  slope.  Thus,  we  can  see  from  this  example 
that  in  the  limit  as  the  number  of  subdivisions  of  the  propagation  interval  tends  to  infinity,  the 
trajectories,  will  for  the  most  part  be  everywhere  continuous;  however,  except  for  the  strait  line 
connecting  source  and  receiver,  not  one  of  the  trajectories  will  be  differentiable  at  every  point  in 
the  propagation  interval.  By  the  very  construction  of  the  trajectories  as  the  limiting  case  of  all 
polygonal  curves  connecting  source  point  to  observation  point,  one  finds  the  set  of  points  where 
the  trajectory  function  to  be  nondifferentiable  to  be  countably  infinite. 


Functional  or  Variational  Differentiation 

There  are  several  definitions  of  functional  or  variational  differentiation  including  the 
GMteaux  variation  and  the  Fr6chet  variation;  however,  we  shall  use  the  following  definition 
(24:22) 

|M  =  lim  +  (R 

Ou[xoj  |Ai|— 0,mai|«u|->0  0u{x)dx 

where  6u,  a  small  variation  or  deviation  in  u,  is  nonzero  only  over  a  small  region  Ax  about  the 
point  xo  and  provided  that  the  limit  exists  for  ail  possible  small  variations  Su  of  u  and 
contractions  of  Ax.  Although  this  definition  appears  to  be  imposing,  it  is  not  appreciably  any 


F-2 


more  difficult  to  use  than  the  definition  function  differentiation  with  which  we  are  more  familiar. 
Furthermore,  an  example  of  functional  differentiation  would  follow  the  first  example  of  function 
differentiation  which  appeared  in  the  previous  section  and  therefore  will  be  omitted. 


Combined  Function  and  Variational  Differentiation 

This  section  contains  some  rules  which  will  be  used  later  in  this  thesis.  These  are  rules 
for  the  variational  differentiation  of  functionals  and  functions  of  functionals.  They  are  given 
without  proof  nor  derivation.  Again,  the  interested  reader  is  urged  to  consult  rigorous  treatment 
elsewhere. 


Rule  1 


S  /  ,r  1  L  r  n 


(F.4) 


where  a  and  b  are  independent  of  u{x). 
Rule  2 


dG[u\  ^  6u{x) 


(F.5) 


where  /  is  a  function  and  G  is  a  functional. 
Rule  3 


6u{x) 


u(y)  =  6{x  -  y) 


(F.6) 


where  u{x)  is  a  function,  6{x  —  y)  is  the  Dirac  delta  distribution  with  argument  (x  —  y). 
Rule  4 

^  ^My)  =  :^r-T-T«(y)  =  -  v) 


Su(x)  dx  dx  Su(x) 

where  S'(x  —  y)  is  the  first  derivative  of  the  Dirac  delta  distribution  with  argument  x-y. 
Rules 

^//(yHy)dy  =  /(x) 


(F.7) 


(F.8) 


where  f(x)  is  a  function. 
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Function  and  Functional  Expansions 


A  function  at  a  fixed  point  may  be  expanded  in  a  power  series  about  another  fixed  point. 
That  is  /(x)  may  be  expanded  about  the  point  y  as 


n=l 


n! 


(F.9) 


where  denotes  the  nth  derivative  of  the  function  f,  Rn  is  the  remainder  after  N  terms.  If  the 
value  of  Rn  tends  to  0  as  N  tends  to  infinity,  the  above  power  series  expansion  is  known  as  a 
Taylor  series.  It  is  implicit  that  the  function  /  and  its  derivatives  are  defined  at  y.  In  operator 
form,  the  Taylor  series  would  appear  formally  as 

/(x)  =i^exp|(x-y)^|/(y)  (F.IO) 


The  equivalence  of  Eq.  F.9  and  Eq.  F.IO  can  be  seen  by  applying  the  usual  expansion  for  exp(x). 


Functional  or  Variational  Taylor  Series 

A  functional  may  be  similarly  expanded  in  a  variational  series.  The  variational  series 
analog  to  a  Taylor  series  is  more  complicated  as  one  would  expect.  Specifically,  a  variational 
Taylor  series  expansion  of  the  functional  F  with  argument  which  is  a  function  u  when  expanded 
about  the  function  v  takes  the  form 


FM  =  FH  + 


(F.ll) 


i  j  -  »(l'))(-(*")  -  »(*"))  + 

^  j  j  J “  v(x'))(u(x")  -  v(x"))(«(x"')  -  vix'")) 

- - dx'dx"dx"'  +  . . .  +  Fyv 

8v{x')8v{x")6v{x"’)  -r  ^  N 
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As  in  the  case  of  the  function  Taylor  series,  if  the  size  of  the  term  Rff  tends  to  zero  as  N  tends  to 
infinity,  the  above  functional  series  is  known  as  a  variational  Taylor  series.  Size,  in  the  case  of 
this  thesis,  will  be  defined  as  the  maximum  of  the  magnitude  of  R^.  For  example,  suppose  we 
have  the  functional  F[x]  =  where  x  is  any  function  from  the  collection  C(o,i].  where  C[o,i)  is 
the  collection  of  all  continuous  functions  defined  on  the  closed  interval  [0, 1].  Clearly,  the 
functions  sin(x)  and  cos(x)  are  members  of  C[o,i].  Therefore,  the  expansion  of  F[sin(x)]  about 
F[cos(x)]  is  given  formally  by 

F[sin(i)]  =  F[cos(i)]+  (F.12) 

(sia(x')  -  cos(x'))(sin(x")  -  + 

^  J I  y  (sinC^O  —  cos(i'))(sin(i")  —  cos{i"))(sin(x''0  —  cos(i”')) 

Scos{x' 

The  application  of  rules  1  through  5  will  result  in 


6^F[cos(x)] 


)6cos(x")^cos(x'") 


dx'dx"dx'"  + ... 


F[sin(x)]  =  F(cos(x)]+  (F.13) 

J  3{s'm(x')  —  cos(x'))cos^(x')^(x  —  x')dx'  + 

1  /  /6(sin(x')  -  cos(x'))(sm(x'')  -  cos(x"))cos(x') 

6(x  —  x')6{x'  —  x")dx'dx"  + 

^  J  j  6(sin(x')  -  cos(x'))(sin(x")  -  cos(x"))(sin(x"')  -  cos(x'")) 

^(x  —  x')S(x'  —  x’')S(x"  —  x"')dx' dx"  dx"'  +  ... 


Simplifying  this  equation  results  in 
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(P.14) 


F[sin(i)]  =  cos^(x)  +  3cos^(i)(sin(x)  —  cos(x))+ 

3cos(x)(sin(x)  —  cos(x))^  +  (sin(x)  —  cos(x))^  =  sin^(x) 


a  tautology,  which  should  not  be  surprising. 
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Appendix  G.  Important  terms 


Action  functional  along  a  path  -  Given  the  trajectory  (also  called  the  path)  of  a  Newtonian 
particle,  the  action  functional  5[x2,  xi]  is  defined  as, 

rh 

5[i2,xi]  =  /  L{x,x,t)dt  (G.l) 

Jit 

where  xi  is  the  initial  location  of  the  Newtonian  particle  at  time  fi,  X2  is  the  final  location 
of  the  Newtonian  particle  at  time  t2,  L{x,  x,  t)  is  the  classical  Lagrangian  of  the  particle. 
The  notion  of  an  Action  functional  can  be  generalized  beyond  that  of  a  Newtonian  particle, 
by  appropriately  generalizing  the  Lagrangian  of  the  particle. 

Feynman  integral  -  An  integral  over  a  function  space.  The  measure  of  the  integral  is  implicitly 
that  of  Richard  Feynman’s  original  construction.  It  is  formally  written  as, 

A'[x2,xi]  =  lini  /  f  •  ■  f  A~'^exp(^S[x2,xi])  Yi  dxm  (G.2) 

<— *07—00  7—00  J—oo  , 

V  II  ’”  =  1 

n— 1  terms 

where  x\  is  the  initial  location  along  the  continuous  trajectory,  5[x2,  xi]  is  the  action 
functional,  A  is  defined  as  =  Planks  constant  divided  by  27r  e  =  is  a  time 

subdivision. 

It  should  be  noted  that  the  Feynman  integral,  unlike  the  sequential  >\lener  integral  with  a 
real  variance  parameter,  does  not  possess  a  countable  measure  and  therefore  is  not 
rigorously  defined.  However,  this  lack  of  rigor  has  not  created  any  difficulties  with  the 
application  of  Feynman  integrals  to  practical  problems,  it  simply  means  that  the 
mathematically  rigorous  rules  governing  integration  in  a  function  space  are  incomplete  at 
present. 

Fraunhofer  approximation  -  A  more  restrictive  approximation  than  the  Fresnel  approximation 
such  that  the  Fresnel  approximation  may  be  approximated  as  a  Fourier  transform  of  the 
optical  disturbance. 

Fresnel  approximation  -  An  approximation  of  the  Huygens-Fresnel  integral  for  small  angles 
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about  the  mean  direction  of  propagation  and  sufficient  distance  from  the  optical  disturbance 
such  that  the  following  approximation  holds. 


(G.3) 


where  \p  —  p',  z\  =  yj{p  —  ff)^  +  z^,  |p  —  =  ^{p-  pf)^,  z  is  the  distance  form  the 

source  of  the  optical  disturbance  to  the  observation  point,  the  integration  region  is 
implicitly  R^. 

Fresnel  zone  -  A  circular  region,  in  a  plane  transverse  to  a  point  source  along  the  mean  direction 
of  propagation.  The  radius  of  the  (first)  Fresnel  zone  is  approximated  as  v/Az,  where  z  is 
the  distance  from  the  point  source  to  the  transverse  plane  and  A  is  the  carrier  wavelength.  It 
originated  from  the  location  of  the  first  destructive  interference  zone  and  from  analysis  of 
weak  propagation  medium  fluctuations.  See  the  figure  below. 


y 


Function  space  -  A  collection  of  functions  which  is  closed  under  (usually)  addition  and 
multiplication  and  has  the  appropriate  identity  elements  and  element  inverses. 

Functional  -  A  mapping  from  a  function  space  into  the  real  numbers.  An  example  of  a 
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functional  is  the  the  mapping, 


X[f\  = 


if  /(x)  is  differentiable  at  0 
otherwise 


(G.4) 


Another  example  is  the  mapping  of  the  area  under  all  real-valued  bounded  continuous 
functions  defined  on  [0, 1]  into  R.  In  this  case,  the  functional  A[f|  could  be  defined  as. 


Mf]  =  /  f{j)dT 

Jo 


(G.5) 


Functional  integral  -  Integration  of  a  functional  over  a  given  function  space.  See  Feynman 
integral  and  >\^iener  integral. 

Green’s  function  -  A  point  source  solution  to  an  associated  operator  equation.  For  example, 
given  L[u(x)]  =  /(x),  the  associated  Green’s  function,  G(x;^)  when  substituted  for  ti(x) 
will  yield:  i[G(x;^)]  =  S{x  —  where,  S{x  —  is  the  dirac  delta  distribution.  For  wave 
propagation,  the  Green’s  function  is  the  spherical  wave  given  by 

Helmholtz  equation  -  A  reduced  form  of  the  wave  equation.  This  equation  takes  the  form, 

( +  k  V(r,  u;))u{r,  a;)  =  /(r)  (G.6) 

where  is  the  three-dimensional  Laplacian  operator,  u(r,  a;)  is  a  (possibly  complex) 
wave  amplitude,  n(r,  u>)  is  the  index  of  refraction,  r  is  a  3-dimensional  position  vector,  k  is 
^  and  A  is  the  carrier  wavelength,  a/  =  ck  is  a  fixed  carrier  angular  frequency,  c  is  the 
speed  of  light,  /(r)  represents  wave  sources  in  the  medium. 

Inhomogeneity  -  A  region  of  a  medium,  such  as  the  atmosphere,  which  exhibits  uniform 
conditions  locally.  These  conditions  may  be  slightly  to  vastly  different  from  the  local 
conditions  only  a  short  distance  away. 

Markov  approximation  -  An  approximation  where  a  stochastic  differential  equation  may  be 
approximated  as  a  markov  process.  This  approximation  enables  the  derivation  of  closed 
differential  equations  for  the  statistical  moments  of  the  independent  variable  of  the  original 
differential  equation. 
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Maricov  process  -  An  independent  increment  random  or  stochastic  process  which  is  conditionally 
dependent  only  on  the  most  recent  interval.  That  is,  the  conditional  distribution  of  all  past 
behavior  intervals  of  the  markov  process  is  dependent  only  on  the  most  recent  interval. 

Maxwell’s  equations  •  A  set  of  electromagnetic  held  equations  which  describes  the  interaction 
and  propagation  of  static  and  dynamic  electric  and  magnetic  fields  as  well  as 
electromagnetic  fields. 

Measure  -  Rule  by  which  differential  contributions  to  an  integral  are  weighted  and  added.  It 
must  be  explicitly  given,  or  in  the  case  of  Riemann  integration  implicitly  given. 

Method  of  smooth  perturbations  -  A  method  used  to  obtain  the  complex  wave  amplitude  in  the 
parabolic  equation.  It  is  restricted  to  weak  intensity  fluctuations  with  the  normalized 
intensity  less  than  0.3.  See  also  Rytov’s  method. 

Parabolic  approximation  of  the  wave  equation  •  Also  known  as  the  parabolic  equation  -  An 
approximation  to  the  Helmholtz  equation  which  results  from  restricting  the  class  of 
admissible  solutions  to  functions  which  vary  slowly  (with  respect  to  the  carrier  frequency) 
in  the  mean  direction  of  propagation.  This  simplification  replaces  the  3-dimensional 
Laplacian  operator  in  the  Helmholtz  equation  with  a  simpler  2-dimensional  Laplacian 
operator  in  the  transverse  plane  and  a  first  derivative  in  the  original  mean  direction  of 
propagation. 

Path  integral  -  A  functional  integral.  This  integral  is  taken  as  the  limit  of  a  multidimensional 
integral  which  ’’samples”  functional  integrand  at  a  finite  number  of  locations.  The  limiting 
process  allows  the  number  of  samples  to  become  unbounded. 

Perturbation  methods  -  A  mathematical  method  of  solving  equations  which  equates  similar 
orders  of  magnitude  in  a  small  perturbation  or  deviation  in  the  original  equation.  The 
equating  of  similar  orders  of  magnitude  in  the  perturbation  term,  allows  the  investigation  of 
zeroth,  first,  second,  and  so  on,  order  changes  in  the  perturbation. 

Polygonal  curve  -  A  piecewise  linear  continuous  curve.  Such  a  curve,  in  general,  is  not 

differentiable  at  each  location  where  two  consecutive  line  segments  intersect  For  example, 
in  Fig.  G  the  derivative  fails  to  exist  for  the  top  most  polygonal  curve  at  locations  a,  b,  c 
and  d.  Let  denote  the  sequence  composed  of  sets  of  points  where  the  polygonal 


curve  is  recursively  subdivided.  Referring  again  to  Fig.  G  let  Si  =  {a,  b,  c,  d}, 

S2  -  {e,f,9,h,i,j,k},Si  =  {l,m,n,o,p,q,r,s,t,u,v,w,x},  aad 
54  =  {A,  B,  C,  D,  E,  F,  G,  H,  /,  J,  K,  L,  M,  N,  O,  P,  Q,  R,  S,  T,  U,  V,  W,  X,Y}.  We  see 
that  by  construction  5i  C  52  C  ^3  C  S4.  We  see  that  in  the  limit  as  n  becomes  unbounded, 
5„  will  in  general  be  a  continuous  curve  which  is  not  differentiable  at  any  rational  number. 
In  other  words,  the  curve  represented  by  lim„_,oo  5„  will  differentiable  everywhere  except 
on  a  set  of  measure  zero.  Such  a  limiting  curve  represents  the  general  differential 
contribution  to  the  Feynman  (path)  integral  just  as  the  differential  contribution 
f{x  +  dx)dx  does  to  the  Riematm  integral  /  f{x)dx. 

Quantum  mechanics  -  Statistical  study  of  kinematics  of  objects  at  the  Fermi  size  (10"*^  cm). 
Objects  which  gready  exceed  this  distance  exhibit  classical  kinematics. 

Rytov’s  approximation  -  Asymptodc  expansion  of  the  complex  amplitude  of  the  opdcal  wave 
that  maintains  only  the  first  term. 

Rytov’s  Method  -  Also  known  as  the  method  of  smooth  perturbations. 

Scale  size  -  Also  known  as:  characteristic  size,  most  commonly  occurring  size,  representative 
size,  expected  size.  In  inhomogeneous  mediums,  the  scale  size  refers  to  the  size  of  the 
random  inhomogeneities. 

Scintillation  -  Observed  intensity  variations  in  an  optical  wave  due  to  random  constructive  and 
destructive  interference  effects. 

Scintillation  Index  -  A  normalization  of  the  intensity  variations.  The  scintillation  index  is 
denoted  by  0^  and  is  defined  as  the  ratio  of  the  intensity  variance  to  the  mean  intensity 
squared. 

Strong  intensity  fluctuations  -  A  condition  in  which  the  scintillation  index  exceeds  1. 

Saturated  Regime  -  Sub-region  of  the  strong  fluctuation  regime  where  the  scintillation  index  is 
asymptotically  decreasing  to  1.  This  region  is  characterized  by  the  almost  totally 
incoherent  intensity  contributions  from  incoming  optical  waves. 

Scalar  wave  equation  -  A  second  order  in  time  second  order  in  space  partial  differential  equation 
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Figure  G.2.  Sequential  evolution  of  a  general  path 


of  the  form. 


(G.7) 


where  T(x,  t)  is  the  wave  amplitude  (possibly  complex),  c  is  a  constant,  x  is  the  space 
variable,  t  is  the  time  variable. 

Stochastic  scalar  wave  equation  -  A  scalar  valued  wave  equation  augmented  with 
non-determiniistic  parameters. 

Stochastic  wave  equation  -  A  vector  valued  wave  equation  augmented  with  non-deterministic 
parameters. 

Taylor  series  -  A  power  series  expansion  of  a  continuously  differentiable  function. 

Weak  intensity  fluctuations  -  No  firm  definition  exhists;  however,  the  generally  accepted 

defination  is  when  the  intensity  variance  normalized  by  the  mean  intensity  squared  is  less 
than  1. 

Wiener  integral  or  sequential  Wiener  integral  -  A  functional  integral  defined  over  Wiener 
measure. 
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